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AN  AXISYMMETRIC 


NEAR  WAKE  ANALYSIS  USING  ROTATIONAL 
CHARACTERISTICS  * 

by 

Mauro  Pierucci 

Polytechnic  Institute  of  Brooklyn 
SUMMARY 

The  near  wake  of  a  cone  in  a  hypersonic  stream  is  analyzed 
by  simultaneously  solving  the  inviscid  region  and  the  viscous  shear  layer  . 

The  inviscid  region  is  solved  by  the  use  of  rotational  axisymmetric 
characteristics.  It  is  assumed  that  viscosity  and  heat  transfer  play  an 
important  role  only  within  a  region  bounded  by  streamlines  which  at  the 
trailing  edge  of  the  cone  are  for  the  most  part  in  the  subsonic  portion 
of  the  boundary  layer.  This  region,  termed  the  shear  layer,  lies 

‘to**' 

between  the  Dividing  Streamline  (or  centerline)  and  'the  Basic  Streamline. 

The  solution  to  the  inviscid  region  is  obtained  by  specifying  conditions 
along  the  characteristic  line  originating  at  the  shoulder  of  the  cone,  and 
by  specifying  the  pressure  distribution  along  a  free  surface  (Basic  Streamline) 
taken  to  be  the  streamline  which  at  the  shoulder  of  the  cone  separates 

^  This  research  was  conducted  under  Contract  Nonr  839(38)  for  PROJECT 
STRATEGIC  TECHNOLOGY  supported  by  the  Advanced  Research  Projects 
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under  contract  DAHC04-69  -C-0077  monitored  by  the  IJ.S.  Army  Research 
Office. 
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Mechanics;  currently,  Research  Specialist,  Electric  Boat  Div.  ,  General 
Dynamics  Corp.  ,  Groton,  Connecticut. 


the  supersonic  from  transonic  and  subsonic  portions  of  the  boundary 
layer.  The  pressure  distribution  along  the  Basic  Streamline  is  iterated 
until  the  mass  flow,  momentum,  and  energy  in  the  shear  layer  are 
consistent  with  the  location  of  the  Dividing  Streamline  and  with  the 
initial  conditions  at  the  edge  of  the  cone. 

Profiles  for  pitot  pressure,  static  pressure  and  stagnation 
enthalpy  are  presented  and  compared  with  experiments  at  different  down¬ 
stream  locations.  The  shape  and  strength  of  both  the  lip  and  recompression 
shock  are  also  shown.  Both  sets  of  results  are  seen  to  be  in  very 
good  agreement  with  the  experimental  results  available. 
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I.  INTRODUCTION 


The  hypersonic  wake  of  both  blunt  and  slender  bodies  has  received 
considerable  attention  within  recent  years;  an  overall  review  of  the 
problem  may  be  found  in  references  1  and  2.  Many  problems  associated 
with  the  far  wake  have  been  analyzed  so  that  the  interest  has  now  shifted 
to  the  solution  of  the  near  wake. 

Chapman^  analyzed  the  problem  of  mixing  of  a  uniform  stream  with 

a  semi-infinite  stagnant  region.  However,  it  was  not  until  later  that  the 

results  of  this  basic  mixing  study  were  used  to  analyze  the  recirculation 

region  and  the  shear  layer  behind  a  blunt-based  body.  In  the  mixing 

process,  a  dividing  streamline  is  obtained  which  separates  the  fluid 

initially  at  rest  from  that  initially  in  motion;  this  streamline,  when  used 

in  conjunction  with  the  actual  body  geometry,  is  assumed  to  divide  the 

recirculation  from  the  external  flow  region  in  the  wake  problem.  Denison 
4 

and  Baum  later  improved  this  analysis  by  solving  the  same  problem  with 
an  initial  (Blasius)  boundary  layer  profile,  which  more  closely  approximates 
the  actual  flow  conditions. 

Few  exact  solutions  have  been  found  to  the  flow  in  the  entire  base 

5 

region;  one  of  these  is  by  Viviand  and  Berger  ,  which  is  valid  for  very 
low  free  stream  Reynolds  numbers.  Their  solution  was  obtained  by  ap¬ 
plying  Oseen's  approximation  to  the  complete  equations  of  motion.  Exact 
solutions  for  laminar  flow  at  higher  Reynolds  numbers  do  not  as  yet  exist. 

Lees  and  Reeve s^  have  attacked  the  near  wake  of  a  blunt  body  by  the 
use  of  the  integral  form  of  the  differential  equations,  as  it  was  done  by  Crocco- 

7 

Lees  ,  and  by  reverse  flow  solutions  to  the  Falkner-Skan  equation.  The  final 
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form  of  the  differential  equations  is  obtained  from  the  x-momentum  and 
from  the  first  moment  of  the  x-momentum  equation.  It  is  assumed  that 
mixing  takes  place  at  constant  pressure  so  that  the  equations  are  sim¬ 
plified  into  two  ordinary  differential  equations  in  two  unkowns  (velocity 
on  centerline  and  displacement  thickness).  Once  the  calculation  is 
carried  to  the  rear  stagnation  point  then  a  new  set  of  ordinary  differential 
equations  is  used  (pressure  is  now  allowed  to  change).  This  new  set  of 
equations  is  now  solved  the  same  way  as  the  previous  ones.  It  turns  out 
that  for  a  given  family  of  solutions  there  will  be  only  one  set  of  values 
which  will  enable  the  calculation  to  go  downstream  (past  the  critical 
point).  For  any  other  values  a  second  stagnation  point  or  zero  pressure 
on  the  centerline  is  obtained.  The  inviscid  flow  field  may  be  assumed  to 
be  governed  by  the  Prandtl-Meyer  equations. 

Due  to  the  vorticity  created  by  the  sudden  expansion  of  the  flow  at 
the  base  of  a  blunt  based  slender  body  in  hypersonic  flow,  the  above 

8 

theory  cannot  be  applied  to  this  class  of  problems.  Reeves  and  Buss 
have  analyzed  this  problem  by  using  the  equations  of  Lees  and  Reeves 
for  the  region  downstream  of  stagnation  point  while  upstream  of  it  the 
Navier-Stokes  equations  are  solved  by  a  double  Taylor  series  expansion 
in  the  stream  function  and  flow  variables  about  the  rear  stagnation  point. 

A  seventh  degree  series  is  used  and  the  coefficients  are  determined  by 
the  symmetry  conditions  along  the  axis,  the  boundary  conditons  and 
temperature  along  the  base  of  the  body  and  the  Navier-Stokes  equations. 
The  outer  inviscid  flow  is  solved  by  the  method  of  streamtubes.  This 
method  may  be  applied  to  two-dimensional  or  axisymmetric  bodies. 
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Rom  and  Rom  and  Victor  have  used  a  modified  form  of  the 
Crocco-Lees  technique  and  with  the  help  of  semi -empirical  results  have 
been  able  to  correlate  experimental  results.  Webb,  Golik,  Vogenitz  and 
Lees**  have  also  extended  the  analysis  of  Crocco  and  Lees  and  have  ob¬ 
tained  results  downstream  of  the  rear  stagnation  point  by  applying  poly¬ 
nomials  in  a  two  moment  (momentum  equation)  plus  one  moment  (energy 
equation)  calculation  and  also  in  a  three  moment  plus  a  two  moment 
system.  This  system  of  resulting  equations  permits  more  degrees  of 
freedom  in  the  choice  of  the  inital  profiles. 

All  the  theories  discussed  above  do  not  lead  to  detailed  solutions  but 

only  give  overall  characteristics  of  the  flow  field  (velocity  on  centerline, 

12 

displacement  thickness,  momentum  thickness  etc).  Weiss  and  Baum 
and  Denison*3  are  the  first  ones  to  have  analyzed  the  flow  field  in  detail. 
Weiss'  analysis  is  limited  from  the  trailing  edge  to  the  rear  stagnation 
point  while  Denison  and  Baum  have  attacked  the  problem  by  starting  at  the 
rear  stagnation  point.  From  trailing  edge  to  rear  stagnation  point,  the 
flow  is  split  up  in  three  regions  (outer  flow,  a  shear  layer,  recirculation 
region).  The  outer  flow  is  solved  by  the  method  of  characteristics.  The 
shear  layer  is  analyzed  by  a  linear  approximation  of  the  boundary  layer 
equations  (Oseen' s  approximation)  while  the  recirculation  region  is  solved 
by  the  inviscid  Navier -Stokes  equations  in  terms  of  vorticity  and  stream 
function.  The  solution  for  the  recirculation  region  is  obtained  by  assuming 
a  temperature  distribution  from  which  a  velocity  distribution  is  arrived  at, 
which  in  turn  is  used  to  solve  the  energy  equation  for  a  new  temperature 
and  pressure  distribution.  Baum  and  Denison  commence  their  analysis 
at  the  rear  stagnation  point  and  integrate  the  equations  by  an  implicit 


3 


A 


» 


difference  scheme.  The  equations  which  are  considered  are  the  continuity 
equation,  the  x-momentum  and  energy  equations  in  the  boundary  layer 
form  and  the  y-momentum  equation  as  applicable  to  an  inviscid  flow. 

The  equations  are  then  integrated  in  the  Von-Mises  coordinates.  How¬ 
ever,  since  the  resulting  equations  for  the  x-derivatives  would  have  a 
singular  point  at  u=a  and  would  be  unstable  for  u<a,  the  transverse 
momentum  equation  for  u<a  is  replaced  by  the  statement  that  p  is  not  a 
function  of  the  stream  function.  (Ihis  replacement  forces  a  physically 
non-existent  saddle  point  on  the  solution).  Now  as  soon  as  any  family  of 
initial  profile  is  picked,  only  one  profile  within  a  given  family  may  permit 
us  to  go  through  the  critical  point.  Any  other  profile  (as  it  happens  in 

Lees-Reeves  theory)  will  give  zero  pressure  or  a  second  stagnation  point 

14 

somewhere  along  the  centerline.  As  explained  by  Weinbaum  Baum  and 
Denison  wrongly  feel  that  if  no  eigenvalue  to  the  particular  family  of  pro¬ 
files  exist,  or  if  two  eigenvalues  exist  then  the  problem  either  is  not  posed 
correctly  or  the  steady  state  solution  as  obtained  from  the  unsteady  equation 
would  have  to  be  analyzed. 

Weinbaum  recently  has  critically  examined  the  differential  equa¬ 
tions  (boundary  layers)  which  have  been  used  by  previous  authors.  He  has 
concluded  that  the  critical  point  obtained  by  most  investigators  is  only  an 
artificial  way  by  which  the  equations  used  (boundary  layer)  manifest  them¬ 
selves  as  not  having  a  dynamic  adjustment  at  the  throat  (i.  e. ,  when 
boundary  layer  equations  are  used  v  at  the  outer  edge  of  the  viscous  region 
cannot  be  arbitrarily  specified,  and  one  has  to  accept  whatever  it  turns 
out  to  be).  Not  only  is  the  critical  point  artificial,  but  its  location  may 
be  varied  at  will  (within  bounds  dictated  by  parameters)  by  suitably  choosing 
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different  positions  for  the  edge  of  the  viscous  region.  The  equations  which 
he  considers  are  the  same  as  the  ones  of  Baum  and  Denison,  except  that 
the  transverse  momentum  equation  is  retained  in  the  subsonic  region. 

The  point  u=a  now  requires  special  care  di£  to  the  fact  that  the  derivatives 
in  the  x  direction  will  be  in  an  indeterminate  form  which  may  be  evaluated 
by  l'Hopital*  s  rule.  With  this  new  set  of  equations  no  eigenvalue  problem 
is  encountered  and  any  arbitrary  set  of  stagnation  point  profiles  will  be 
able  to  pass  downstream.  The  correct  solution  will  then  be  obtained  when 

the  ambient  pressure  is  recovered  at  the  end  of  strong  interaction  region. 

15 

The  Rudman-Rutin  equations  are  similar  to  the  ones  used  by  Weinbaum 
with  the  exception  that  the  x  derivative  of  the  pressure  term  has  been 
neglected.  This  minor  difference  causes  a  major  breakthrough  in  the 
solution,  because  the  singularity  at  u=a  has  now  disappeared  and  the 
numerical  technique  used  may  be  simplified  considerably. 

The  flow  field  described  above  is  not  amenable  to  a  single  solution 
unless  the  complete  Navier -Stokes  equations  are  utilized.  A  solution 
can  also  be  obtained  by  splitting  the  problem  into  the  following  four 
distinct  flow  regimes:  1)  leading  to  trailing  edge  of  body,  2)  expansion 
of  fluid  at  the  trailing  edge,  3)  trailing  edge  to  rear  stagnation  point 
4)  rear  stagnation  point  to  downstream  infinity.  For  the  sake  of  sim¬ 
plicity  it  has  to  be  assumed  that  no  interactions  between  the  different 
regions  take  place. 

If  the  flow  is  assumed  to  be  laminar  from  the  leading  to  trailing 
edge  of  the  cone,  then  the  solution  can  be  obtained  for  either  of  the  two 
conditions.  For  large  Reynolds  numbers  the  inviscid  flow  field  is  solved 
by  either  the  method  of  characteristics  or  by  solving  the  inviscid  conical 
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flow  equations.  Once  the  inviscid  field  is  known  then  the  thin  viscous 
layer  around  the  body  can  be  solved  by  the  usual  boundary  layer  techniques. 

I 

For  low  Reynolds  numbers  the  problem  is  more  complicated  because  no 
distinct  viscous  layer,  shock  wave  and  inviscid  regions  exist  and  the 
viscous  layer  (even  to  a  first  approximation)  cannot  be  neglected  with 
respect  to  the  inviscid  region. 

The  expansion  of  the  rotational  flow  at  the  trailing  edge  is  a  prob¬ 
lem  which  is  now  being  studied.  This  problem  is  complicated  by  the  up¬ 
stream  influence  of  the  base  pressure  through  the  subsonic  portion  of  the 

boundary  layer.  This  problem  has  been  investigated  by  Weinbaum^  for 

17  18 

incompressible  flow,  Baum  and  Weiss  and  Nelson  for  supersonic 
high  Reynolds  number  flows.  In  all  the  cases  mentioned,  no  interaction 
between  the  boundary  layer  and  inviscid  flow  field  is  assumed. 

For  the  analysis  to  proceed  from  the  trailing  edge  to  the  rear 
stagnation  point,  the  recirculation  region,  shear  layer  and  the  outer  in¬ 
viscid  flow  would  have  to  be  solved  independently  remembering  that  the 
boundary  conditons  connect  the  three  solutions  together  (obviously  for 
low  densities,  this  procedure  cannot  be  adapted  because  of  the  interaction 

problems  involved).  The  recirculation  region  and  the  shear  layer  can 

19  20 

best  be  analyzed  by  the  methods  developed  by  either  Weiss  or  Moretti  . 
Weiss'  method  while  not  as  detailed  as  the  approach  used  by  Moretti, 
has  the  advantage  of  being  soluble  within  a  short  period  of  time.  In  both 
cases  the  inviscid  outer  flow  is  solved  by  the  method  of  characteristics. 
Moretti  has  numerically  solved  a  modified  form  (viscosity  is  retained 
while  for  simplicity  heat  transfer  is  neglected)  of  the  unsteady  Navier- 
Stokes  equations.  Due  to  the  hyperbolic  nature  of  the  equations  a  Lax- 
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Wendroff  technique  is  used  to  obtain  the  solution.  The  steady  state 
solution  is  then  assumed  to  be  the  asymptotic  time  limit  of  an  unsteady 
flow  field.  By  using  these  equations,  both  the  recirculation  and  shear 
layer  region  can  be  solved.  With  slight  modifications  the  equations  could 
also  be  used  for  low  density  flows.  However,  the  usefulness  of  this 
method  is  offset  by  the  enormous  time  required  to  obtain  a  solution  (on 
an  IBM  7094  the  time  for  an  accurate  calculation  would  be  of  the  order  of 
several  hours.) 

A  semi-empirical  approach  which  can  be  worked  out  with  the  aid  of 

the  method  of  characteristics  and  the  equations  used  by  Weinbaum  and 
21 

Garvine  or  the  ones  originally  analyzed  by  Rudman  and  Rubin  will 

now  be  outlined.  The  main  tool  to  be  used  in  this  analysis  will  be  the 

method  of  characteristics.  Application  of  the  method  of  characteristics 

to  the  near  wake  was  first  suggested  by  M.  H.  Bloom  in  a  presentation 

$ 

at  an  I,  D.  A.  Conference  in  1963.  Calculations  showing  the  importance 
of  radial  pressure  gradients  and  the  thickening  of  the  shear  region 
(after  the  expansion  of  the  surface  boundary  layer)  immediately  downstream 

of  the  shoulder  of  an  axisymmetric  body  were  also  shown  by  Bloom 

22 

and  Vaglio-Laurin“  .  The  first  published  results  of  this  method  applied 

12  29  23 

to  the  near  wake  problem  are  by  Weiss,  Weinbaum  '  and  Weiss 
24 

and  Weinbaum  .  In  reference  12,  the  base  region  of  the  flow  over  a 
wedge  is  treated  and  an  approximate  solution  is  obtained  by  matching  the 
free  shear  layer,  recirculation,  and  inviscid  flow  regions.  The  assumptions 
of  both  Chapman,  and  Denison  and  Baum  that  the  stagnant  (recirculating) 
region  is  semi-infinite  is  no  longer  necessary  and  thus  the  effect  of  finite 


*  l*he  proceedings  of  these  meetings  are  unpublished. 
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base  diameter  is  obtained. 

In  reference  29  the  variation  of  the  entropy  within  the  boundary 
layer  was  studied  and  in  reference  23,  preliminary  analyses  for  characteristic 
calculations  are  initiated.  It  is  also  shown  that  for  high  "inviscid"  Mach 
numbers  (Mg  >  8) ,  less  than  half  the  total  free  stream  expansion  occurs 
in  the  centered  expansion  at  the  corner  .  The  remainder  of  the  expansion 
is  produced  by  the  reflected  'vaves.  To  show  this,  the  problems  of 
interaction  of  a  slip  stream  with  a  weak  expansion  wave  and  also  the 
interaction  of  a  shear  layer  with  a  wake  expansion  fan  are  solved.  In 
reference  24,  preliminary  calculations  from  the  characteristics  program  are 
presented. 

The  present  paper  presents  a  method  which  combines  the  rotational, 
axisymmetric  characteristics  with  a  viscous  inner  region,  to  determine 
near  wake  profiles.  Imbedded  shocks  are  considered  in  the  characteristics 
solution.  The  surface  boundary  layer  profiles  at  the  separation  point  provide 
the  initial  conditions  for  the  characteristics  program  in  the  supersonic 
region,  while  the  subsonic  part  of  the  boundary  layer  is  taken  into  account 
by  dividing  this  portion  of  the  boundary  layer  into  strips  and  considering 
each  strip  to  be  governed  by  the  one-dimensional  flow  equations  including 
viscosity  and  thermal  conductivity.  The  heat  transfer  and  shear  acting  on 
each  streamtube  are  computed  from  the  average  values  of  temperature  and 
velocity  in  each  strip.  Details  of  the  recirculating  flow  region  are  not 
considered  in  this  analysis.  The  present  analysis  is  useful  to  evaluate 
the  flow  field  with  reasonable  accuracy  to  a  few  base  diameters  downstream 
of  the  body.  At  this  location,  the  profiles  could  then  be  used  as  initial 
data  to  the  available  far  wake  analyses. 
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II.  INVISCID  ROTATIONAL  FLOW  FIELD 


In  order  to  analyze  the  near  wake  (cf.  Fig.  la)  short  of  using 
the  complete  set  of  fluid  mechanical  equations  available  (Navier- 
Stokes),  certain  simplifying  assumptions  are  made.  For  the  present 
analysis  they  are: 

1.  A  steady  state  solution  is  assumed. 

2.  No  interaction  from  the  subsonic  part  of  the  boundary  layer 
(this  alters  the  initial  profiles  and  can  be  readily  included  if  a 
more  accurate  determination  of  this  effect  is  known.  ) 

3.  Expansion  of  the  first  streamline  (Basic  Streamline)  takes 
place  by  means  of  a  Prandtl-Meyer  fan  (A-B  in  Fig.  lb). 

4.  Basic  Streamline  (originally  this  streamline  has  a  Mach 
number  equal  to  M^  at  the  trailing  edge  of  cone)  is  a  free 
streamline  and  its  shape  is  determined  by  assuming  a  specified 
pressure  distribution  along  it. 

With  the  above  assumptions,  once  all  the  variables  are  specified 
along  a  first  family  characteristic  line  emanating  from  the  point  of 
boundary  layer  separation  on  the  cone,  where  the  Mach  number  in  the 
boundary  layer  M  =  M^  >  1.  0,  the  inviscid  flow  field  may  be  analyzed  by 
the  method  of  rotational  axisymmetric  characteristics  including  imbedded 
shocks.  The  manner  in  which  the  pressure  distribution  along  the  basic 
streamline  is  specified  will  be  described  later. 

As  the  flow  near  the  shoulder  expands,  it  will  separate  from  the 
base  of  the  cone  and  a  Basic  Streamline  (B.  S)  is  formed  which  will  separate 
the  inviscid  outer  flow  from  the  viscous  layer  which  is  obtained  by  expanding 
the  subsonic  portion  of  the  boundary  layer.  As  the  external  streamlines 
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progress  downstream,  they  will  reach  a  point  where  their  velocity  will  be 
in  the  direction  of  the  axis.  After  this  point  is  reached,  the  streamlines 
will  again  start  curving  outward  and  the  compression  waves  formed  by  their 
divergence  will  coalesce  into  a  trailing  shock  which  may  be  weak  or  strong 
depending  on  the  cone  angle  and  flow  conditions. 

The  details  of  the  shoulder  expansion  region  are  better  shown  in 
Fig.  (lb).  Note  that  a  lip  shock  may  be  formed  by  the  concave  curvature 
of  the  B.  S.  It  is  also  noteworthy  to  mention  that  in  order  to  solve  for  the 
entire  supersonic  region  one  would  have  to  reflect  the  incoming  expansion 
waves  from  the  sonic  surface.  However,  this  is  not  done  due  to  the 
complexity  involved  in  the  determination  of  the  sonic  surface  which  is 
imbedded  in  the  viscous  region,  therefore,  a  streamline  with  an  initial  Mach 
number  M1  (Basic  Streamline)  is  used  as  the  reflection  surface. 

The  problem  is  mathematically  well-defined  once  all  the  conditions 
along  the  initial  characteristic  line,  the  point  about  which  the  Prandtl-Meyer 
expansion  occurs,  and  the  subsequent  basic  streamline,  are  specified. 

The  analysis  may  then  be  subdivided  into  four  separate  unit  problems: 

1.  Evaluation  of  an  interior  point  (3)  once  values  at  (1) 
and  (2),  are  known  (see  Fig.  2). 

2.  Reflection  of  a  second  family  characteristic  line  from  a  pressure 
surface  whose  pressure  is  a  given  function  of  x.  (Fig.  3) 

3.  Evaluation  of  conditions  at  a  point  which  is  obtained  as  a  result 

of  two  characteristic  lines  of  the  same  family  intersecting  (Fig.  4). 

4.  Extension  of  a  shock  wave  once  conditions  at  a  point  behind  the 
shock  are  known  (Fig.  5)  (i.  e. ,  how  to  obtain  C  once  A  and  B 
are  known). 

To  solve  the  first  unit  problem  the  following  equations  are  available: 
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along, 


^  =  tan  (  0  +  U) 


(first  family  characteristic) 


cotH  JsinSsinji  ,  cosHsin8^  dS  JV  1  cosM  dH  „ 

“T“  dV  '  de  ■  YcosTeTn)  dX  +  "y"Kcos  (0+m)  Sfl  dX  ‘  T*  w¥I9+m)  3TT  dX  =  0 


along, 


dY 

XX  =  tan  (6  -Ui) 


(second  family  characteristic) 


cotji  ,,,  ,  ,A  Jsin0sin(i  ,v  coslisin2U  SS  ,  .  1  cosu  dH  _ 

T"  dV  +  de  '  Ycos(d-n)  dx  ‘  yRcos(S-h)  TO  dx+T8  cos(9 -|i)  3TFT  dX  -  ® 


along, 


tan  0 


(streamline) 


and  also 


S  =  const. 


H  =  const. 


/  M3  + 


All  the  variables  may  be  nondimensionalized  as  follows: 


,  H 
h  =  -rj—  8 

Hoo 


V  X  „  Y 

VL  x=  TT  r-  X 


the  differential  equations  are  then  reduced  to  a  form  amenable  to  solution 
by  a  computer,  and  they  are 


*  =  ^ 


> +  ji°ba 

b„  -  a. 


+ 


=  ha  - 


N  -  gl 

(s3  -  )  ma  AXa 


nx  Axx  +  m3  Axg 
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v3  =  jAt  vx  +  Aa  va  -  (9X  -03  )  +  Bj  A*!  +  d  a  AXg  [c^n^  [3%  -  c^  Ax^ 

+  A$1  ■'  <5«> 

o3  =  ea  -a,  (va  -  v.)  +  d ,  ax,  +  (.,-■! )  flXa  -  (b.  -h, 

All  coefficients  are  defined  in  the  List  of  Symbols. 

The  equations  for  the  reflection  of  an  expansion  line  from  a  pressure 
surface  are  derived  by  assuming  that  the  B.  S.  passes  through  a  point  (1) 

(see  Fig.  3  )  and  conditions  at  a  point  off  the  streamline  (2)  are  known;  the 
continuation  of  the  streamline  is  desired  and  this  is  done  by  locating  point  (3). 

A  second  family  characteristic  line  from  2  to  3  and  a  streamline  from 
1  to  3  are  used.  The  equations  available  are  then: 


»3  =  8i  (6c) 

h3=hi  y-_l_  (6d) 

M®3  =  7H  jf1  +  **0  (  _1  (6e) 

v3  = 

03  =  0a  +  Aa  (va  -  v3^j^  )  +d8&«8  -  ca  (  s3  -  sa)  -  ^  ^  )  (6f) 

where  all  the  functions  are  evaluated  in  the  above  order. 
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In  the  case  of  two  intersecting  characteristic  lines  of  the  same 
family  (AC  and  BD  'a  Fig.  4)  a  shock  is  assumed  to  form  at  the  inter¬ 
section  point  E.  By  using  the  values  at  E_(on  line  AC)  and  E+(on  line  BD), 
the  strength  of  the  shock  is  found.  Since  it  turns  out  that  the  shock  is  a 
very  weak  one  (i.  e. ,  -^2L  <  2-3%),  the  line  BD  is  extended  by  using  the 
first  family  characteristic  line  in  front  of  AC  as  the  line  from  which  a 
second  family  is  eminated. 
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III.  SHEAR  LAYER 


* 


* 


The  portion  of  the  boundary  layer  whose  Mach  number  is  subsonic 
before  the  expansion  at  the  trailing  edge  cannot  be  analyzed  by  the  method  of 
characteristics  since  the  shear  and  heat  conduction  effects  may  be  large  in  this 
region;  in  addition,  a  portion  of  this  layer  remains  subsonic  even  after 
expansion  at  the  corner.  It  is  therefore  analyzed  by  conserving  mass 
momentum,  and  energy  in  individual  stream  tubes  wherein  the  flow  is  assumed 
to  be  one  dimensional,  including  the  effects  of  shear  and  heat  conduction 
between  different  stream  tubes.  This  procedure  also  determines  the  pressure 
distribution  along  the  outer  most  stream  tube  which  is  adjacent  to  the 
characteristic  field.  By  matching  the  two  regions,  therefore,  the  analysis 
can  proceed  downstream  in  a  consistent  manner. 

The  subsonic  part  of  the  boundary  layer  is  divided  into  "n"  strips  and  it  is 

assumedthat  each  strip  expands  inviscidly  and  adiabatically  from  pg  to  p^ 

as  if  it  were  one  dimensional.  As  soon  as  this  corner  expansion  is  completed, 

each  streamtube  is  followed  by  using  the  one  dimensional  equations  with  shear 

25 

and  heat  transfer  as  given,  for  example,  by  Shapiro  (see  Fig.  6).  The 
energy,  momentum,  and  continuity  equations  in  nondimensional  form  are 


dh. 

i 

-sr 


=  2tt 


M 

t  _ 

M.  A.p. 

i  li 


hi  yi  ‘  ViViJ 


(7) 


dM.3  , 

nr-  =_2m.®  (i  +  -ij— 


dp.  .  dh.  XT. 

’  I  i  i 


Mi3  }  [yp-kfo  (  “ET  }  $.w  "cEr  + 


1  1 


^i  -yi-i 


(8) 


d K. 

-£T  =  *1 


(1-M.3)  dp. 


YPiM.3  dx 


+  (1 


(1  +  (Y  - 1)  M.3)  c, 

.  .  dh.  T. 

+  Y^i  M.a)  *  _L  + - 1 

l  h.  dx  v  •  -  v  •  * 

l  1 l  yi-l 


(9) 


where 


Qi  = 


K 


3T 

■Sy 


(T.  -  T.  .) 
i  l-r 

fri-'i-i1 


Va  c 
°~Z~  Cf. 


av  AVi 
M  Ty  ■  TTY” 


The  unknowns  in  the  above  equations  are  A^,  hL,  h^.  The  pressure 
p.(x)  is  assumed  to  be  known  and  is  equal  to  the  value  as  given  by  the 
inviscid  characteristic  program. 

Zero  heat  transfer  is  assumed  along  the  basic  streamline  and  a 
mild  stagnation  temperature  variation  is  assumed  in  the  recirculation 
region.  Internal  diffusion  takes  place  by  virtue  of  heat  transfer  and 
shear  across  the  strips.  Since,  for  a  given  pressure  distribution  p(x), 
the  basic  streamline  has  been  obtained  from  the  characteristic  program, 
this  line  is  used  as  a  reference  line  from  which  the  radial  dimension  ox 
the  n  strips  is  measured  when  the  "correct"  pressure  distribution  is 
used  along  the  basic  streamline.  The  boundary  line  of  the  inner  most 
pressure  distribution  used  to  determine  the  characteristic  field  and  the 
location  of  the  basic  streamline  is  therefore  iterated  on  to  obtain  this 
condition  before  the  analysis  proceeds  to  the  next  streamwise  calculation 
of  the  matching  flow  fields.  In  this  manner,  the  dividing  streamline  and 
the  rear  stagnation  region  can  be  obtained  once  the  base  pressure  p^ 
is  specified. 
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IV.  RESULTS 


* 


Few  detailed  experimental  profiles  of  flow  parameters  in  the  near 
wake  of  blunt  based  cones  are  available  in  the  literature.  Most  experimental 
papers  dealing  with  this  subject  present  results  in  terms  of  variables  which 
do  not  describe  the  details  of  the  flow  field.  For  example,  there  is  a 
significant  amount  of  information  available  on  base  pressure,  heat  transfer 
to  base  or  rear  stagnation  point  location  while  detailed  measurements  of 
the  flow  field  are  not  presented. 

Two  recent  papers  which  present  local  profiles  of  pressure,  temperature, 

26 

etc. ,  at  various  downstream  stations  are  by  Schmidt  and  Cresci  and 
27 

Bauer  .  The  free  stream  Mach  number,  Reynolds  numbers  and  other 
pertinent  test  conditions  are  presented  below  for  the  two  experiments. 
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The  above  conditions  were  used  as  inputs  for  the  characteristics  program. 

In  both  cases,  the  experimented,  values  of  the  base  pressure  were  used 
in  conjunction  with  the  axisymmetric,  rotat  mal  program.  Since  this 
information  can  be  obtained  from  empirical  correlations  (cf.  Ref.  28,  for 
example)  for  different  flow  conditions,  no  generality  is  lost  by  this  assumption. 
Once  the  base  pressure  is  known,  the  pressure  distribution  from  the 
base  to  the  rear  stagnation  point  can  be  determined  if  the  maximum  Mach 
number  (minimum  pressure)  is  specified.  As  explained  in  the  previous 
sections,  the  location  of  the  rear  stagnation  region  is  obtained  by  matching 
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the  viscous  shear  layer  with  the  inviscid  characteristic  program.  The 
viscous  shear  layer  from  the  cone  base  to  the  rear  stagnation  region 
was  found  to  be  governed  principally  by  the  effect  of  heat  transfer  from 
the  recirculation  region  while  downstream  of  the  stagnation  region,  internal 
shear  produced  the  largest  effect.  The  temperature  of  the  layer  adjacent 
to  the  shear  layer  was  assumed  to  be  either  (i)  a  constant  or  (ii)  vary  from 
the  cone  surface  temperature  at  x  =  0  to  the  recirculation  region  temperature 
at  the  rear  stagnation  point.  As  seen  from  Fig.  (9-a)  the  effect  due  to  this 
variation  appears  to  be  negligible.  For  the  Mach  eight  conditions,  profiles 
at  different  x  stations  are  obtained  and  values  of  pitot  pressures,  static 
pressure  and  stagnation  enthalpy  are  plotted  and  compared  with  the  experi¬ 
mental  results  in  Figures  (7)  through  (  9).  The  stagnation  enthalpy  profiles 
are  seen  to  agree  very  well,  especially  in  the  region  close  to  the  rear 
stagnation  point.  The  accuracy  of  the  analysis  decreases  in  the  downstream 
direction  as  the  region  dominated  by  diffusive  effects  grows  into  the  flow 
field  computed  by  characteristics,  thereby  invalidating  the  basic  assumption 
of  an  inviscid  outer  flow.  The  pitot  pressure  profiles  are  seen  to  be  in 
good  agreement  up  to  x/D  =3.25.  In  contrast  to  the  other  two  sets  of 
profiles,  the  static  pressure  profiles  are  less  accurate  close  to  the  rear 
stagnation  point.  This  is  believed  to  be  due  to  two  effects.  First,  it  is 
much  more  difficult  to  accurately  measure  static  pressure  in  the  recirculation 
and  stagnation  region  due  to  probe  interference,  and  second,  the  theory  is 
able  to  determine  the  local  pressure  distribution  at  every  point  downstream 
of  rear  stagnation  point  in  a  self-consistent  manner,  while  in  the  rear  stagnation 
region  the  uniqueness  of  the  static  pressure  distribution  is  not  guaranteed. 


It  is  also  seen  that  at  each  x  station,  the  location  of  the  lip  and 
recompression  shock  can  be  predicted  with  good  accuracy.  In  Fig.  (10) 
one  sees  the  flow  field  for  this  case.  The  light  lines  indicate  the  various 
first  family  characteristics  emanating  from  the  basic  streamline  (for  clarity, 
second  family  lines  are  omitted),  while  the  two  darker  lines  show  the 
ocation  and  shape  of  the  lip  and  recompression  shock. 

In  Fig.  (11)  pitot  pressure  profiles  for  the  Mach  3.  0  case  are  presented. 
Due  to  the  few  results  published  by  Bauer,  the  only  variables  which  were 
compared  were  the  pitot  pressure  profiles;  again  comparison  between  theory 
and  experiment  is  reasonably  good  and  the  location  of  the  lip  shock  is 
predicted  within  experimental  accuracy. 

One  may  note  that  the  Crocco-Lees  type  of  singular  behavior  doesn't 
appear  in  this  analysis  since  the  critical  region  is  not  analyzed  in  detail. 


V.  CONCLUDING  REMARKS 


The  present  analysis  of  the  near  wake  is  able  to  predict  (with 
reasonable  accuracy)  flow  conditions  and  shock  shapes  at  different  locations 
without  the  necessity  of  analyzing  the  recirculation  region  in  detail.  The 
advantage  of  this  is  immediately  evident  in  that  the  problem  of  the  recir¬ 
culation  region  is  formidable  due  to  the  complexity  of  the  differential 
equations  which  govern  the  flow  and  the  specification  of  boundary  conditions 
along  an  undetermined  boundary.  It  may  therefore  be  inferred  that  unless 
one  is  interested  in  the  recirculation  region  per  se,  the  detailed  solution 
to  this  region  will  not  play  an  extremely  important  role  in  the  downstream 
flow.  The  complete  solution  of  the  recirculation  region  has  been  replaced 
by  the  specification  of  several  conditions  in  this  region:  1)  the  base  pres¬ 
sure  (ref.  23),  2)  an  average  temperature  of  the  recirculation  region 

(ref.  24),  3)  maximum  Mach  number  (minimum  pressure)  along  the  center - 

line  (  ref.  25).  Since  these  parameters  empirically  have  been  correlated 
under  different  flow  conditions  in  the  referenced  papers,  these  conditions 
can  be  readily  obtained. 

The  following  conclusions  may  be  drawn:  1)  stagnation  enthalpy 
profiles  are  relatively  insensitive  to  the  shape  of  the  initial  profiles  and 
to  the  heat  transfer  from  the  recirculation  region,  2)  while  this  is  also 
true  for  the  inviscid  portion  of  the  stagnation  and  static  pressure  profiles, 
the  shear  layer  is  quite  sensitive  to  these  conditions. 

As  a  result,  further  refinements  and/or  extensions  to  the  present 
work  should  deal  with  1)  detailed  analysis  of  the  initial  expansion  of  a 
compressible  shear  layer,  2)  analysis  of  the  rear  stagnation  region,  and 


19 


* 


* 


3)  a  better  representation  of  the  heat  transfer  and  shear  between  the  shear 
layer  and  the  recirculation  region. 

The  freedom  in  inputs  of  the  present  theory  will  allow  for  calculations 
of  near  wake  of  more  general  shaped  bodies  (spheres  or  blunt  bodies).  The 
only  modification  involved  would  be  the  alteration  of  inputs  along  the  character¬ 
istic  line  emanating  from  the  separation  point. 
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FIG.  (4)  INTERSECTION  OF  TWO  FIRST  FAMILY  CHARACTERISTIC  LINES 
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FIG.  (5)  PROCEDURE  USED  TO  EXTEND  SHAPE  OF  SHOCK 
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FIG.  (7)  PITOT  PRESSURE 
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FIG.  (8)  STATIC  PRESSURE  PROFILE,  M^'8.0 
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APPENDIX  A 


DESCRIPTION  OF 'CHARACTERISTIC  PROGRAM 

The  program  is  capable  of  analyzing  the  trailing  edge  expansion 
and  the  near  wake  for  any  given  cone  angle,  Mach  number,  and  free  stream 
stagnation  conditions.  It  consists  of  13  different  sections;  eight  of  these 
are  function  subroutines  for  the  different  coefficients,  the  remaining  five 
are  MAIN,  CINPUT  {  Calculate  INPUT  ),  CHAR  (  CHARacteristic  ),  CSHOCK, 
(  Cone  SHOCK  ),  DIVST  (  Dividing  STreamline  ),  subroutines.  A  flowchart 
for  each  of  these  is  provided  by  Figs.  A-l  to  A- 5. 

At  the  start,  the  MAIN  program  directs  the  computer  to  the  CINPUT 

subroutine.  The  function  of  this  subroutine  is  to  evaluate  and  store  in  the 

memory  of  the  computer  the  initial  characteristic  line  and  the  conditions 

along  a  streamline  which  originates  from  a  point  where  the  Mach  number  is 

equal  to  M^  ,  The  conditions  along  this  streamline  change  by  means  of  a 

Prandtl-Meyer  expansion  procedure.  In  order  to  evaluate  the  initial  first 

family  characteristic  line,  the  "viscous"  part  of  the  line  is  solved  and 

matched  with  the  "inviscid"  region.  The  inviscid  or  potential  characteristic 

22 

line  is  read  in  to  the  computer  from  Sims'  tables  .  This  line  is  then 
shifted  a  little  upstream  and/or  downstream  until  the  two  lines  match  to 
the  desired  degree  of  accuracy.  To  evaluate  the  viscous  part  of  the 
characteristic  line,  the  only  inputs  which  are  necessary  are  various  values 
of  Mach  number  varying  from  M^  to  Mg.  Once  the  Mach  number  is  known, 
the  velocity  u/ug  can  be  found  (all  the  values  at  the  edge  of  the  boundary 
layer  are  assumed  to  be  equivalent  to  the  inviscid  values  on  the  cone  and  are 
obtained  from  Sims'  tables).  Once  u/ug  is  known,  corresponding  values  of 
x  and  y,  may  be  found  by  using  Blasius  solutions  for  a  cone  or  by  assuming 


45 


any  desired  boundary  layer  profile.  More  features  of  this  and  subsequent 
subroutines  may  be  obtained  from  the  flow  diagram,  and  the  program. 

From  the  CINPUT,  the  MAIN  program  calls  the  CHAR  subroutine  which  is 
used  to  evaluate  conditions  at  a  third  point  once  the  conditions  at  two  other 
points  not  on  the  same  characteristic  are  known  (i.  e. ,  this  is  just  a  rotational 
axisymmetric  characteristic  program  where  the  rotationality  comes  from 
the  boundary  layer  profile). 

Once  all  the  points  along  a  given  expansion  line  are  known,  the  shape 
of  the  shock  originating  from  the  tip  of  the  cone  may  be  found;  this  procedure 
is  carried  out  in  the  CSHOCK  subroutine. 

When  the  expansion  fan  is  completed  the  program  goes  on  to  the  DIVST 
subprogram.  This  subroutine  evaluates  conditions  on  the  dividing  stream¬ 
line  by  using  previous  first  family  characteristic  lines  that  the  computer  has 
evaluated.  This  having  been  done,  the  computer  goes  back  to  the  CHAR 
program  to  evaluate  the  next  point  on  this  new  first  family  line. 

After  the  expansion  fan  is  completed,  it  will  be  seen  that  characteristic 
lines  of  the  same  family  will  tend  to  cross  each  other.  Whenever  this 
happens  the  location  of  the  intersection  is  found  and  checked  to  see  whether 
the  assumed  shock  is  strong  or  weak.  If  the  shock  is  strong  enough,  formation 
of  an  imbedded  shock  is  postulated,  the  necessary  values  upstream  and  down¬ 
stream  are  evaluated  and  the  program  goes  on  to  the  next  characteristic 
line  in  the  flow.  If  it  is  weak  shock  (as  determined  by  the  As/ s  criteria) 
evaluation  of  the  characteristic  line  is  continued  with  a  new  reference  line 
used  to  evaluate  the  rest  of  the  points  on  this  line. 

After  every  characteristic  line  is  completed,  at  given  values  of  x, 

M,  p/Pqq.  H,  p  /pt  ,  and  y  are  evaluated  so  that  profiles  of  these  valuesvs. 
y  may  be  plotted  at  different  downstream  locations. 
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MAIN  PROGRAM 


FIG.  (A-l)  INVISCID  PROGRAM 
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SUBROUTINE  CHAR 


FIG.  (A- 3)  INVISCID  PROGRAM 
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SUBROUTINE  CSHOCK 


FIG.  (A- 4)  INVISCID  PROGRAM 
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SUBROUTINE  DIV  ST 
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FIG.  (A- 5)  INVISCID  PROGRAM 
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COMPUTER  PROGRAM  FOR  CHARACTERISTIC  CALCULATION 

COMMON  V(  5,70)  ,  THETA  (  5, 70  )  ,  AM  (  5 , 70  )  ,  S  (  "> ,  7C  )  ,  X  (  5 , 7C  )  ,  Y  ( 5 , 7C  ) ,  NL  ( 5  ) 
1  ,EPI (5, 70), AH (5, 73) , PM (5, 70) ,L0C(5) ,XP(5,70) ,XX(5, 150), 

1  GAM , RS , AM  I , N , I , J , M , J AM , AL P , T B E T , NP , K A , M AXN , MN , MC , I M , N PT , L F SH , 
lLL,TL,Pi,EMl,YT,IDEL,IRED,KK,AL,A4,B4,C4,C4,E4,F4,IR,Ch,SP,PSI,TSI 
DIMENSION  UGUE  ( 4  1  )  *  ETA ( 41  )  ,XZ(  10) 

DIMENSION  P2X( 13) , P 2 P ( 1 0 ) 

KEWiNDl 

REWIND2 

READ(5,370)GAM,RS,AMI,ALP,TBET 
READ (5,380)  IN 
READ (5,380)  INF 
READ ( 5  »  3fcO )  PI, EMI 
READ ( 5 , 340 ) IRED, IDEL 
READ ( 5 , 520 )  YT 
READ ( 5 , 520  >  XT 
RE AD ( 5 , 3B0 ) J  AM 
READ  (5, 350  LL  ,  TL 
RE AD ( 5 , 380 )  LN 
READ ( 5 , 490 )  ( XZ ( K ) , K=1 , LN ) 

READ(5,380)LPSH 
READ! 5, 380 )  KK 
IF(IN.EQ.l)  GO  TO  10 
RE AD ( 5 , 380 ) N 
READ ( 5 , 380 ) MAXN 

READ(5,390)(V(1,J),THETA(1,J),AM(1,J),S(1,J),X(1,J),Y(1,J),EPI(1,J 
1  ) , AH (1,J) *  J  =  1  *  N ) 

GO  TO  20 

10  CALL  CINPUT 
MAXN=N+KK- 1 
20  INN=N 

I MAX=MAXN 
DO  30  J-2  ,  N 

3'«  PM  ( 1 ,  J  )  *  (  Y  (  l*J)-Y(l,J-i)  )/(X(i*J)-X(l,J-l)  ) 

G2=(GAM-1. )/2. 

UN1=1.+G2»AMI**2 
PX=GAM/ (GAM-1.  ) 

40  KA  =  1 

IR  =  1 
Ns  INN 
MAXN= I  MAX 
NRED=5 
MC=0 
NP*  1 
M*0 

RE  AD  (5, 480  A4 , 84 , C* , D4 , E4 , F4 

DL=AL/. 0174532925 

TOC=TS I *Gw 

SRT= 12 • *SR 

PSID=PSI/144. 

WR I TE ( 6 ,4 33)  AMI,TSI,PSIC 
IF  C JAM. EO. 1 )  WRITE(6,410)  DL 
IF(JAM.EO.J)  WR  I T  E ( 6 ,420 )  DL 
WR I TE ( 6 , 430  )  SRT  ,  TCC , A4 , EM  1 
1  =  1 
MN=0 
I M=KK- 1 
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5'1  WRITE(6»**5J) 

WK  IT:  (  6,4fc '•  )  IR 

WRI  TE  (6,A7'1)  (  V(  I  ,  J)  ,THCTA(  I  ,J)  ,  AM(  I  ,  J)  ,S(  I  ,  J)  ,X(  I  ,  J)  ,Y(  I  ,  J  )  ,  AF 
1  (  I  ,  J  > , J=i,N) 

W  R  I  T  0  (  6 , 4  4  •'* )  E  P  I  (  I  »  N  ) 

WRITE(6»51-'  ) 

RO  50  LX=1,LN 

IF(XZ(LX).LT.X(  I  tm  i,0  TO  80 
00  6C  L  J  =  1 »  5  0 

I  F  (  X  (  I.LJI.5Q.  j.  )  GC  TO  8" 

I F  (  XZ  C  LX )  .  L  T . X (  I  ,LJ)  )  GC  TC  70 
60  CONTINUE 

7u  FL  =  ( X (  I , L  J ) -XZ ( L  X )  )/(X(  I,LJ)-X(  I,LJ-1)) 

FH= ( XZ ( LX ) -X (  I »  L J- 1  )  )/(X(  I  ,  L J  )  -  X  (  I.LJ-l)) 

T V  =  THE  T  A (  I ,  LJ-1 ) *FL+ THETA (  I , LJ) *FF 
AMV=AM (  I iLJ-1 ) • FL+ AM {  I,LJ)*FH 
TEN!P  =  ( 1 .  +  ARV  * • 2/ 5  .  ) 

SV=S(  I ,LJ-1 ) *  FL ♦ S (  I »  L J ) • Fh 

VV=V(  l »  L J- 1 ) • FL  +  V (  I,LJ)*FH 

YV  =  Y (  I »  L J- 1 ) *FL ♦ Y (  I »  L  J ) * FH 

HV=  AH  (  I  t  LJ-1  )  +FL+ AH (  I,LJ)*Fh 

PV=(UNl/( l.+G2*ARV**2)  )**PX*EXP(-SV)*FV**PX 

SIG  =  PV«TEKP*VV*YV*COS (TV J/SORT (FV  ) 

PT2=(  (GAM+1.  )/2.*ARV**2/UM)**PX*(  (  G AR ♦  1 .  ) /2 .  /  ( G  AR*  AR V*  *  2-G 2  )  ) 

1  **(1. /(GAM-1. ) )*PV 

WRITE (6, 500)  XZ(LX),IR,YV,SIG,ARV,VV,SV,PV,PT2,HV 
WRITE(2)XZ(LX),I«,YV,SIG,ARV,VV,SV,PV,PT2,HV 
80  CONTINUE 

IF(X(I,l).oT.XT)  CALL  EXIT 
NIU  I»=N 
1  =  1  +  1 
I  R  =  I  R  +  1 

I F ( I . GT . 3 )  GC  TC  290 
9  J  =  1 

N  =  N  ♦  1 

IF  (Nl.GT.RAXN)  GC  TC  150 
10..  IF(IR.GT.KK)  GC  TC  16) 

IF ( INF. EC. 1 )  GC  TC  28  j 
IFII.GT.2)  GC  TO  11 J 
CM  =  GAM- 1 . 

GP=G AR+  1 . 

TF  =  THETA(1 , 1  ) 

AMF=AM(  1,1) 

ARK  *S0RT(2./GM*(  CM»*Z*GR/?.  +  l.  )#(  ( PI/ A*  )  ••(GM/GAM  )  )-l.  )  ) 

CNUL  =  S0 U (GP/GR ) *  AT  AN ( S&R T ( GR/ GP* ( ARK **2- 1 .  )))-ATAN(SCRT(ARK*«2. 
1  -1.  )  ) 

CNUF  =  SCR  HGP/GF  )  *ATAN ( SCRT (CR/GP* ( AR (  1  ,1)*«2-1.)))-ATAMSCRT 
1  (  A  f*  (  1  ,1)** 2.-1.)  I 

TKK  =CNCF+TFETA(1,1 J-CNUL 

XKK=KK-1 
SAT=SC«T(OP/GR  ) 
llu  X(I,1)=X(1,1) 

Y(I,1)=Y(1,1) 

S( I, 1)=S(  1,  1  ) 

AH ( I ,  1  )  = AF ( 1 ,  1  ) 

EPI ( 1 , 1  )=CDI ( 1  ,  1  ) 
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AMU  ,  1  )  =  AMF 
X IK= IR-1 

THCT  A ( I » 1 )  =  TF  + (TKX-TF )*XIK/XKK 
CNU=CNUF  +  TF-THt:TA{  T  t  1  ) 

DC  120  K I =1 » 100 

CNW-SAT*ATAN(SCRT < ( AM (  I  , 1 ) **2-1 . ) *GM/GP ) ) -AT AN ( SCRT { AW (  I  ,  1)«*2 
1  -1.1) 

AMU  *  1  )  =  AM  {  I  , 1  )  +  2.*(CNL-CNl*  ) 

IF(ABS(CNL-CMn/).LE..0uC01)  GC  TC  13C 
12.  CONTINUE 

130  VII  V1)«SQRT(AK(I  , 1 )**2/ ( 2 . /GM+AM (I  ,1)«*2)) 

140  CALL  CHAR 
GO  TO  5c 
150  N  =  MAXN 
MN  =  1 

GO  TC  10^ 

160  IF(NREO.EG.O)  GO  TC  170 
CALL  OIV  ST 
IF(N.EQ.l)  GC  TO  4C 
I F ( N. EQ . 2  )  GC  TO  270 
GO  TO  140 

170  I F ( IRED.LE.N)  GO  TC  180 

IREDMREC-IDEL 
GO  TO  170 
180  ILP=IR-3 
REWIND1 
L  =  5 
LB=0 

DO  250  ILC0P=1 » ILP 
READ(1)IW,NL(5),L0C(5) 

NU I =NU ( 5 ) 

READ(l)(V(5»J)»THETA(5»J)fAM(5tJ)»S(5tJ)tX(5tJ)  ,  Y ( 5 , J  ) , EP 1 1  5 , J  ) , 
1  AH(5fJ)fPM(5tJ)«XP(5tJ)fJsl«NUl) 

READ(1)(XX(5,IX)»IX=1,IW) 

BACKSPACE  1 
BACKSPACE l 
BACKSPACE  1 
190  K=1 

DO  200  J=l,  IRECt  ICEL 

V  (L,K)  =  V  (  L  » J  ) 

THETA(L.K)*THETA(Lf J) 

AM  ( L ♦ K ) = AM  ( L  t  J  ) 

S  (L,K)  =  S  ( L  t  J  ) 

X  ( L  »K  )  =  X  (L»J) 

V  ( L  t  K  )  =  Y  IL»J) 

AH  ( L  » K  )  = AH  (L»J) 

EPI  ( L  *K 1 =EPI  (L»J) 

20o  K=K+1 

INEX=IRED+1 
N  =  N  U  (  L  ) 

I F ( INEX.GT.N)  CC  TC  22: 

DO  210  J=INEX,N 

V  (L»K)=V  ( L  *  J  ) 

THETA(L»K  )  =  THETA(L»J  ) 

AM  ( L  »  K  )  =  AM  ( L  t  J  ) 

S  ( L  » K  )  *S  (LtJ) 
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X  (  L  »  K  )  =  X  (  L  •  J  ) 

Y  (  L  *  K  )  =  Y  (  L  •  J  ) 

AH  (  L  »  K  )  =  AH  (  L  •  J  ) 

FPI  ( L  *  K )  =  E  P  I  (L,J) 

l1  ;  K  =  K  +  1 
i?'  NU  (  L  >  =  K  —  1 

N  =  NU  ( L ) 

OG  2  3u  J  =  2  *  N 

23  J  PM!L,J)=!Y(L,J)-Y(L,J-1))/(X(L,J)-X(L,J-1)) 

OC  ?4j  IX=1,IW 
240  XX(L» IX)=1C.E+10 

IF (L.EC.S  )  GC  TC  26  3 
250  CONTINUE 
LB=Lb+  I 
L  =  LB 

IFIL.LT.3)  GC  TC  190 

NRED=1 

MAXN=N 

1=3 

GO  TC  16^ 

260  WRITE! 1 ) IW.NO! 5) ,LCC(6) 

WRITE{l)lVl5,J),TF5TA(5,J),AM5,J),S(5,J),X(5,J),Y(5yJ),EFI(5,J), 
1  AH(5,J),PM(5*J)*XP!5»J),J=1,N  ) 

WRITE(l)  ( XX ( 5  » I  X ) *  IX=1,IW) 

GC  TC  25C 

270  IF(MC.NE.C)  CALL  SHOCK 
IF(MPT.fcO.l)  CALL  CHAR 
GO  TC  50 

2B0  RE AD (5, 390)  V ( 1 , 1 ) , T HE T A (I , 1 ) , AH ( I , 1 ) , S ( I , 1 ) , X ( I , 1 ) , Y ( I f 1  )  , E P  I  (  I , 

1  1 ) i AH ( I , 1  ) 

GC  TO  140 
290  I  W=  l  R-3 

NU I =NU (  1  ) 

WRITE(l)  I W , NU  I  *  LCC (  1 ) 

WRITE{l)(Vll,J),THFTA(l,J),AMl,J),Sll,J),X(l,J),Yll,J),EFI(l,J), 
1  AH ( 1  *  J )  ,PM(1,J),XP!1,J),J=1,NUI  ) 

WRITE(l) ( X  X (  1 ,  IX),  I  X  =  1 ,  I W  ) 

CO  310  1=1,2 
N=NU ( I + 1 ) 

DO  3CJ  J  =  1 , N 

Y  ( I , J)=  V  (  I  +  l, J) 

THETA! I, J)=  THE  T  A ! I  +  1 , J ) 

AM  l I , J)=  AM  ! I  +  1 ,  J  ) 

S  (  I  ,  J  )  =  S  !  I  +  1 ,  J  ) 

X  !I,J)=  X  ( I+l, J) 

Y  l I , J ) =  Y  (  1  +  1, J) 

EPI  l I , J )=  e  P I  ( I+l, J) 

AH  ( I , J)=  AH  (  I+l,  J) 

PM  ( I , J  )  =  PM  l  I  +  l, J) 

30C  XP ! I , J ) =XP (  I  ♦  1 ,  J  ) 

LOCI  I ) =LGC l  I  +  l  ) 

NU! I ) =NU {  I  +  l  ) 

I  W  =  I  W  +  1 

DC  310  I  X  =  1 ,  IW 
3 1 C  XX! I , I  X ) =  X  X  l  1  +  i  ,  IX  ) 

DO  320  I  X  =  1 ,  I R 


55 


320 


* 


XX ( 3  , 1 X ) =0 . 

1*3 

N=NU ( 3  ) 

00  333  J  =  1  *  N 
V  (  I  ( J  )  *  c  • 

THETA ( I v  J )=G  • 

AM ( I , J ) =0  . 

S!I,J)=0. 

XIIiJI-0. 

Y( I , J ) =G. 

AH ( I »  J ) =0  • 

XP{ ! ,J)=0. 

PM ( I , J ) *G • 

333  EP I ( I »  J ) =0. 

GO  TO  90 

340  FORMAT (215) 

350  FORMAT ( 1 3 ,  E 1 3. 6 ) 

36'J  FORMAT  (  2E 18  •  8 ) 

370  FORMAT ( 5E1 5 • 6 ) 

380  FORMAT (12) 

390  F0RMAT!4E18.8/4E18.8) 

400  F0RMAT(/////»1X,25HFREE  STREAM  MACH  NUMBER  = , F18 . 8 , // , IX , 36HFREE 

1  STREAM  STAGNATION  TEMPERATURE  = , F 18 . 8 , 3X , 15HCEGREES  RANK INE , //IX , 

1  33HFREc  STREAM  STAGNATION  PRESSURE  =  ,  F 1 8. 8 , 3X » 3FPS I ,  /) 

410  FORMAT! IX »2CFHALF  ANGLE  OF  CCNE  =, F 18 . 8 , 3X , 7hCEGREE S ,  /) 

42w  FORMAT! IX, 21HHALF  ANGLE  OF  WEDGE  * , F 18 . 8 , 3X , 7FCEGREE S ,  /) 

430  FORMAT (IX,24HRA0IUS  OF  BASE  CF  CCNE  * , F 18. 8 , 3X, 6HINCHES » / / , IX , 2 1H 
1TEMPERATURE  CF  CONE  * , F 18 . 8 , 3X , 15FDEGREES  RANK  INE , //IX , 30HFRES SLRE 
1  BASE  /  PRESSURE  INF  *,F18.8,//,1X, 4 9H INITIAL  MACH  NUMBER  CF  TRAIL 
1ING  EOGE  STREAMLINE  =,F18.8) 

443  F0RMAT!100X,10HEPSILCN  *  ,E18.8) 

45C  FORMAT!///, 124H  V  THETA  M 

IS  X  Y 

1H  ) 

460  FORMAT ( /// 8H  ROW  ,12///) 

470  FORMAT ( 7E1 8 • 8 ) 

480  FORMAT ( 6E 13 . 6 ) 

49 0  FORMAT! 1QF7. 3) 

500  FORMAT (2X,2HX=,F7.3,1X,  I2,4X,8E13.4) 

510  FORMAT!///, I IX , 3HR0W , 13X , 1HY , i IX , 3FS IG , 1 IX , 1HM , 1 2X, 2FVV , 1 IX , 

1  2HSV,10X»4HP/PI,9X,7HPT2/PTI»7X»2FHV) 

520  “FORMAT ( E18.8 )  ”  ... 

END 


SUBROUTINE  CINPUT 

COWON  V(5,70)  .THETA  (5, 70),  AM  5,70i  ,S(5,70)  ,  X  <  5 , 7C  )  ,Y(5, 7C  )  *M  (5) 
I  ,EPI  (5,  ?<*  )»AH(5,7':  ),PM(b,70)  ,LGC(5)  ,XP{5,70)  ,XX(  5,150, 

1  CAM,  PS  ,  AM  ,M,  I  ,  J  ,f' ,  JAM,  AL P  ,  T e E T , N P , K A , V AXN , MN , MC ,  IMNPT.LPSF, 
1LL,TL,P1,EV1,VT,ICEL,IR^U,KK,AL,AA,B4,CA,C4,E4,F4,IR,CK,SR,cSI,TSI 
DIMENSION  UCLE (41 ) ,ETA( 41  ) 

PE  AO (5, 19 D  ( r  T  A ( I ) ,UGUF  (  I  ),  1  =  1,41) 

READ!  5,2Cj)EM  ,GV>,TST  ,PSI  ,R,CPP,CI ,  AL,SR,TCT,RGR,AMIS 
READ!  5,24C)AP 
R  E  AC ( 5  »  2  20 )  IY 
GM  =  GAV-1  . 

GP=GAM+1. 

XL=SR/SI\ ( AL ) 

PI=1.+GM/2.*AVI**2 

PE=l.+GN/2.*FME**2 

TE=TSI*TCT/PI 

RE=RCR/PI*»(1./GM)«PSI/(P»TSI) 

UE  =  AMS*SCRT(GM/GP*2.*CPP*TSI) 

EU*2.27E-8*TE**1.5/( 19e.fc+TE ) 

REL=RE*UE*XL/EU 

XL=XL/SR 

AJA= JAM 

RR=2.*XL/SQRT( (2. «AJA+1 . )*REL) 

UEU  =  cME/AM»SO«T  (PI/PE  ) 

HRITl(6,?30)  RtL,Rl:,6U,XL,Ut,UEC 

READ(5,220)  NC 

REAO( 5, 24C  )  ( AH( 1, J ) , J=1 ,NC  ) 

READ ( 5 , 220 )  NL 
J  1  =  N  C  +  1 
J2=NL+\C 

REAC(5,210)(X(1,J»,Y(1,J),TH£TA(1,J),S(1,J),AM1,J),AM1,J),J  =  J1, 

1  J2) 

DG  1C  J=J1,J? 

EPI (  1  ,  J  )  =  : . 

V  ( 1 ,  J)=SQRT(GM/GP)*APM1,J) 
n  AM( 1 , J) =  S«RT ( (2./GM*V( 1 , J)**2) /( 1.-VC1, J )**2> ) 

READ(5»24G)  E  P  I  (  1 , J  2 ) 

VE  =  SQRT ( 1 ./ ( l.  +  2./GK/EME**2)  > 

HRITE(6,299)  Gw,Vt 
DO  80  J  =  1  , NC 
EPI  (  1 , J ) =  0 . 

V( 1, J)=SCRT(l./(l.  +  2./GM/AMl»J)**2) ) 

VEV=( VE/V( 1, J) ) **2 

UU=  (  ( 1  .“GM  +  SCKT  (  (  l.-GW)**2+4.*Gfc*VEV)  )/(2.*VEV) 

WRITE(6,2H0)  UL 
IF(IY.EC.l)  GG  TC  IB" 

IF1UU.LT. UCUEtl ) )  GC  TO  30 
IF(UU.GT.UGLE(41) )  GC  TO  30 
DO  20  1=1,41 

IFlUU.tC.lClEl I ) )  GO  TO  4C 
IFIUU.LT.U0UE1 I ) )  GO  TO  5C 
20  CONTINUE 
3CT  HR  I T  b ( 6 , 2  5C )  UL 
CALL  EXIT 
40  T  A  =  P  T  A  (  I  ) 

GO  TO  6 C 
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50  TA*ETA(  1-1)4  (  E  TA  (  I  ) -ETA  (  I- 1 1  J/IUCUE  (  d-UCLE  (  1-1)  )  *  ( Ul-UCLE  ( I  - 1 )  ) 

60  FX»1.-EXP(-TA) 

fcXX=l.-EXP(-2.*TA) 

AY  »RR*(GH*PE*TA  +  PE*  d . -GW ) *( TA-C I *EX )- 1  PE- 1 .  )*( TA-2.*C I*EX  + 

1  CI**2/2.*EXX) ) 

Y(l, J)=AY*CCS(AL)4l.O 
70  AH(l,J)«Gk+(l.-GK)*UU 

THETAd,J)*AL 

Sd,J)*-2.*GAM/GM*(AL0G(AMd, J  ) /EME  ) -AL0G  UU  )  )«S(1. Jl) 

UX*U ( AM ( 1 , J ) ) 

AFF*AF 

AF*COTAN ( THETAd»J)+UX) 

IF ( J.EO. 1 )  GC  TO  170 

X(1,J)*X(  l,J-l)  +  (AF4AFF)/2.*(Yd,J)-Yd,J-l)  ) 

80  CONTINUE 

90  IF  (XU,ND).GT.X(  1,  J2))  GO  TO  160 

00  100  Js J1 , J2 

IF(X(1,ND).LT.X( 1, J) )  GO  TO  110 
100  CONTINUE 

110  ZK*Y(l,J-l)4(Y(l,J)-YCl,J-in/(Xd,J)-X(l,J-i)  )*(Xd,NC)-Xd,J-l)  ) 
BE«Y(1,ND)  -  ZK 
WRI TE ( 6 , 270 )  BE 
00  120  J*J1,J2 
Y(1,J)*Y(1,J)*(1«40E ) 

120  XU,J)=XU,J)M1.4BE) 

IF(ABS(BE).GT..C3CC005)  GC  TC  90 
00  130  J*J1,J2 

IF(Xd,J)»GT*Xd,NC) )  GO  TO  140 
130  CONTINUE 

140  L*J1-1 

DO  150  K»J,J2 
L«L4l 

X(1,L)«X(1,K) 

Y ( 1 ,  L ) *Y ( 1 , K ) 

V(1,L)*V(1,K) 

S ( 1 , L ) aS ( 1 ,K ) 

AH (1,1)  *AH ( 1 ,K ) 

AM(l,L)«AP(l,K) 

EPI ( 1*  L )*EPI ( 1 ,K ) 

THETA(1*L)«THETA(1,K) 

1 50  CONTINUE 

N»ND4J2-J4l 

IF (  IY.EC.O)  MR  I TE ( 6, 30C )  NC 
IF(IY.EQ.l)  WR I TE ( 6 , 310 )  NC 
RETURN 

160  WRITE ( 6,260 )  X ( 1 , NC )  ,  X ( 1 , J 1 ) , X ( 1 , J2 ) 

CALL  EXIT 

170  Xd,l)*-AY*SIN(AL) 

I F ( IY.EC.1)  X(1,1)=-(UU**AP*,066*SR*EPE**.824)/REL**.116 
GO  TO  80 

180  Yd,  J)*(LL**AP*.066*SR*EPE**.824)/REL«*.  116*C0TAN(AL)4l.C 
GO  TO  70 

190  F0RMAT(F6.4,E13.6) 

200  FORMAT  (4E18.8) 

210  FORMAT ( 3E 18 • 8 ) 

220  FORMAT ( 12) 
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2  30  FORMAT!  LH1,  1X,4HREL=  , G 1 5 . 8 , 5H R HOE = , E 1 5 . 8 , 4HMIE = , E 1 5 . E , 2HL=  ,E15.8, 
1  3HUE=,£15.8,6HLE/UI=,£15.8) 

24*  FORMAT  (E18.fi) 

250  FORMAT (IX*  9HL  GVPR  UE  ,  E  18 . 8 , 1 2HCLT  GF  RANGE) 

260  FORMAT ( 1X»15HLAST  DARK  PC  I N T , E 1 8 . 8 , 24HN0T  BETWEEN  LIGHT  POINTS, 

1  2618.8) 

27u  FORMAT  !2X*3FBE  =  »E18..8) 

2b0  FORMAT  ( 2 X , 3HUU  = , E 1 8 . b  ) 

*90  FORMAT!/, IX , 3HGW  = , E 1 8 . 8 , 3HVE  =  , E 1 8 . 8  ,  / / ) 

300  FORMAT!/////, IX,  26HLAVINAR  BOUNDARY  LAYER  - ,I5,16H  POINTS  ARE 

1  USED) 

310  FORMAT!///// , IX, 28HTURBULENT  EOUNCARY  LAYER  - ,I5,16H  POINTS  ARE 

1  USED) 

END 
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SUBROUTINE  CHAR 

COMMON  VI  5,70  .THETA  (5,70)  ,  AM  5, 70)  tSI  5,70)  ,X(5, 70)  ,Y  (5,70  ,NL(5) 
1  ,EPI (5,7^), AH (5,70) , PM (5, 70) ,L0C(5) ,XP(5,70)  ,  XX  I  5,150 t 
1  G AM »  RS  *  AM  I ,N,I , J,M, J AM , ALP , TBE T , NP , KA , M AXN  ,MN , MC , I M , N PT , L P SH , 

1  LL,TL»P1,EM1,YT, I  DEL , I  RED ,KK , AL , A4 , B4 ,C4  ,D4 , E4 , F4, 1 R 
T (XX,YY)=COS(XX)/CCS (YY+XX) 

Z(XX,YY)=COS(XX) /CCSIYY-XX) 

WRITE  (6, 520  (  XX  (I  - 1 , 1  2 )  ,  I  2  =  1 ,  I R  » 

IFINPT.EC.il  GO  TO  140 
G0=0. 

LIM=N-1*M 

IF(MC.GT.O)  LIM=LIM-1 
INTEGERP 

DO  330  J  =  2  » L  I  M 

10  IF(GC.EC.l..ANC.XP(I»J-l).EQ.O.)  GO  TC  120 
IF (GO.NE. 1. )  GO  TO  130 
IQ=  I  I 

IFIII.LT. (IR-2) )  GC  TO  40 
IF(  (  IR-m.EC.2)  I  X*1 
IF ( (  IR-I I ) .EC. 1 )  I X  =  2 
20  IQ=IQ-1 

I F ( IO.EG.O)  GO  TO  400 

IFIXXI Itl I l.GT.XX(IXflO) )  GO  TO  20 

I X*  I Q 

I  I  =  I  X 

I F ( I  I.LT. (IR-2) )  GC  TO  50 
I F ( (IR-II J.EC.2)  I X= 1 
IF(UR-II).EO.l)  IX  =  2 
30  P*J+1 

IFIII.LT. KK)  P=P+ I I-KK 
GO  TO  140 
40  I WHE - 1 

GO  TO  60 
50  I WHE  =  2 

60  LBACK=0 

70  BACKSPACE l 

BACKSPACE  1 
BACK  SPACE  1 

READ(l)  I S » NU ( 5  ) . LOC ( 5 ) 

IF(II.EO.IS)  GO  TO  80 
BACKSPACE  1 
LBACK=LBACK+1 
GO  TO  7C 
80  NU I  =NU ( 5 ) 

READ!  1MV(5,K)  ,THETA(5,K)  , AM  I  5 , K ) , S < 5 ,K )  ,X(5,K  )  ,  Y  ( 5 ,  K  )  ,  EP I  I  5  ,K  )  < 

1  AH(5fK),PM(5,K),XP(5»Kl,K*l,NUn 
READ (1) (XXI5.IX) ,IX«1* IS) 

I F ( L BACK . EQ .0 )  GO  TO  100 
DO  90LBA= 1 . LBACK 
READ(l) ITT  »NU(4) » L  0  C (4) 

NU I =NU 1 4 ) 

READ(l) (V(4,K) ,THETA(4,K) , AM (4, K ) , S I  4 , K ) ,X(4,K) tY(4,K) ,EPI(4,K  )* 

^  1  AHC4»K) »  PM ( 4  » K ) »XP{4»K)»K31,NUI ) 

90  READ! 1) I X X ( 4 » I  X ) , I  X  =  1 » ITT) 

100  I X=5 

IF! IWHE.EC. 1)  GO  TC  20 
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GO  TC  30 
113  IXM  -1 
1 I  =  I R- 1 
GO  TO  170 
12.'  P  =  P+1 

GO  TO  140 
130  I  X  = I  -1 
1 1  =  I R- 1 
P=J+NP-l 
14u  L=0 

WRITEI6,51C)  IR,J,II,P 
150  L=L+1 

160  IFINPT.EG.l)  GC  TC  110 
170  IF(l.EQ.i)  GC  TO  450 
180  HG=HX-GX 

XI  I ,J)  =  (Y(I, J-l)-YIIX,P)+HX*X< IX,P)-GX*X< I,J-1) )/HG 

Y(  I ,  J)  =  (GX*HX*(X( IX,P)-X( I,J-l i J+HXfrY (I , J-1)-GX#Y< IX,P) )/FG 

OELX=X( I *  J ) -X I I X  *  P  ) 

DELZ=X I I , J)-XI I, J-i) 

XN=DELZ*ANX 
XM=DELX*SMAX 
SS=S(IX  ,P)-SII,J-i) 

DELH=AH I IX  ,P)-AHI I , J-l  ) 

XNM=SS/IXN+XM) 

DX=OW*DELX 
S  3  =  S ( I »  J  ) 

SII,J)=S(  I  » J-l )+XNM*XN 

AH ( I , J ) -AH  I I ,  J-l ) ♦ I (DELH/IXN+XM) )*XN) 

S1=V ( I ,J) 

IF(AH(I»J).LT.O.)  AH<  I  ,  JIM  AMI-1,  J)+AH(  I-1,J  +  1)  )/2. 

190  V(I,J)S(AX*V(I , J-l I ♦ AW* VI IX,Py+THETA( IX,P)-THETA(I,J-1)+BX«DELZ 

1  +DX+SS*ICX*DELX»SMAX-DELZ*ANX*CW)/(XN+XM)+DELH/<2.*IXN+XM) )*<TX* 

1  DELZ/I AH (I, J-l )*V< I, J-l )**2)~ZW*CELX/(AH(IX,P)*V( IX,P)**2) ) )/IAX 
1  *SQRT ( A  H ( I , J ) /AH ( I ,J-1 ) )  +  AW*SQRT  f  A  H ( I ,J)/AM IX, P)  )  ) 

S2  =  THE  TA (I , J ) 

THETA!  I  ,  JIMHETAI  IX,P)*AW*(V<IX,P)-SCRT<  AH(  I  ,  J)/AH(  I  X,  P  )  )«Vt  I,J)  ) 
1  +DX+CX*XM*XNM-ZW*DEIX*DEIH/I2.*VI IX,P)«*2*AH( I X ,P ) * l XN+XM ) ) 
IFIS2.NE.0..ANC.L.EQ.1)  GO  TC  250 
IF(V(I,J).LT.O..OR.V(I,J)«GT.l.)  GO  TC  390 
AMI  I, J)=SQRT(2./(GAM-1. ) ) * SQR T C VII, J 1**2/ (1.- VII , J)**2) ) 

I F ( ABS ( ( V  <I,J)-S1)/V  ( I , J ) ) •IT*  TL )  GC  TC  240 
200  EPI(I,J)=0. 

PM  II , JIMYt I , J ) -Y  C I, J-l) )/IX( I,J)-X( I, J-l) ) 

IF(IR.LE.KK)  GO  TO  220 
210  IF(L.GT.LL)  GO  TO  420 

220  I F I  AM  I  I , J  )  • LT* 1 • I  AM ( I , J )  =  ( AM  I  1-1 , J ) ♦ AM ( I- 1 , J  + 1 )  ) /2  . 

IFIL.EQ.l  )  GO  TO  460 
230  UX=UI AMI  I  ,J)  ) 

AX=.5*(AZ+A(tX,V( I , J ) ) ) 

AW=.5*IAY+A(LX»V( I ,J) )  ) 

DW=.5*(DY  +  0( JAM, UXV THETA  I  I, J),Y(  I,J)  ) ) 

BX=.5*(BZ+B( JAM,UX,THETA(I,J),Y(I,J)) ) 

CW=.5*(CV+C(UX,GAM,RS) ) 

CX  =  .5*ICZ+C(UX,GAM,RS)  ) 

'  ZW=.5*IZZ+Z(UX,THETA(I,J) )  ) 

TX  =  «5*(TZ  +  T(CX, THETA! I,J)  )  ) 
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GX=.5*{GZ  +  TAN! THETA (I  ,J  )+UX)  ) 

HX  =  .5#(HZ  +  TAN(THETA(ItJ)-UX)  ) 

IF( IR.LE.KK.ANC.L.GT.LLS  GO  TO  310 
GO  TO  150 

240  IF ( V ( I , J ) *GE  •  1  •  )  GC  TO  210 

IF (ABS { (THETA ( I , J ) -S2 ) / THET A (I , J ) ) . GE  .TL  .OR . ABS ( ( S ( I , J ) -S3 ) /S ( I , J 
1  )).GE.TL)  GC  TO  210 

IF(Y(I,J).LT.O..CR.Y(I,J).LT.TBET)  GC  TO  400 
250  PM(I,J)=(Y(I,J)-Y(I,J-1))/(X(I,J)-X(I,J-1)) 

IF(L.EQ.l)  GC  TO  260 
IF(IR.LE.KK)  GO  TO  310 
IF(PM(I,J).GT.PM(IX  ,P  ))  GC  TO  260 
GO  TO  310 

260  YPM  =  ( PM ( I • J  )«Y ( I X  ,P-1)  -  PMIIX  ,P)»Y(I,J-1)  4  PM(I,J)« 

1  PM  ( I X  , P ) •( X ( I » J-l )  -  X(IX  ,P-1)  ) )/  (PM(  I  ,J)  -PM(IX  ,P)) 

XP(I,J)=X(I, J-1)4((X(I,J)-X(I,J-1) )/(Y(I,J)-Y(I, J-l) ) ) * ( YPM-Y I  I  ( J 
1  -1)1 

IF(L.EQ.l)  GC  TC  270 

I F ( X  P ( I t  J ) • L  T • X (  IX  tP))  GO  TC  270 

XP(I«J)*0. 

GO  TO  3 1C 

270  WRITE(6»480)IRfJ»YPM 
L0C(I)*64L2 
IFCL2.LE.0)  L2=l 
MC=2 

IF(XP(I,J).LT.X(I, J-l) )  GO  TC  370 
WRITE ( 6»490  )  IR,J,XP(I,J) 

XX(I,II ) *XP ( I • J  ) 

YA*Y(I ,J-l)-YIIX,P-l) 

YB*Y( IX, P)-Y(IXtP-l) 

YC=YPM-Y( IX, P-1) 

XA=X( IX,P)-X( IX, P-1) 

XB=X ( I ,J-1)-X( IX, P-1) 

xc=xp(i,j)-y  nx, p-i) 

DEN=XA*YA-Yt?*XB 

DFXT  =  ( (THETA  1  I  X , P ) -THETA ( I  X , P-1 ) ) *YA-YB* (THETA ( I , J-l )-THETA(  IX, 

1  P-1 ) J ) /DEN 

DFXH= ( (AH(IX,P)-AH(IX,P-1))*YA-YB*(AH(I,J-1)-AH(IX,P-1) ))/CEN 
DFXV= ( ( V  (IX,P)-V  (IX, P-1) )*YA-YB*(V  l  I ,  J- 1 ) -V  ( I X , P-1 ) )  ) /DEN 
DFXS= ( ( S  (IX,P)-S  (  IX, P-1) )*YA-YB*(S  ( I  ,  J-l ) -S  ( I X, P-1 ) ) ) /DEN 
DFYT=(XA« ( THETA ( I , J-l) -THETA ( IX, P-1 ) )-XB*( THETA ( I X , P )-THETA 
1  (  IX, P-1) )  )/DEN 

DFYH=(XA*(AH( I »  J- 1 )-AH ( IX, P-1 ) )-XB*(AH(IX,P)-AH( IX, P-1) ) )/CEN 
DFYV  =  (XA*(V  (I,J-1)-V  (IX, P-1  )  J-XBMV  ( I  X ,  P ) -V  (  IX,  P-1 ))  ) /DEN 
DFYS= ( X  A  * ( S  ( I , J- 1 ) -S  (IX, P-1) )-XB# (S  ( I X , P )-S  ( IX, P-1 ) )  )/CEN 
THETA(IX,P)=THETA( IX, P-1 ) 4QFXT*XC4DFYT*YC 
V( IX,P)=V(IX,P-1)4DFXV*XC4DFYV»YC 
S( IX,P)=S(IX,P-1 )4DFXS*XC4DFYS*YC 
AH(IX,P)*AH(IX,P-1  )4DFXH»XC4CFYH#YC 

AM ( I  X  , P)=SCRT(2. /(GAM-1. ) ) *  SORT ( V ( IX  , P  )  #*2/ ( 1 .-V ( I  X  , P  )  «*2 ) ) 
YA=Y ( I,J)— Y( I,J— 1) 

YB=Y( I X , P-1 ) -Y (  I, J-l) 

YC  =  YPM-Y (  I  ,  J-l  ) 

XA»X( IX.P-1 )-X( I , J-l ) 

XB*X ( I , J ) -X ( I, J-l) 

XC  =  XP( I , J )  —  X ( I ,J-1  ) 
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DEN=XA#YA-YB*XB 

OF  XT  =  ( ( THETA(  I  X.P-l ) -THETA ( I  ,J-1) ) *YA- ( THET A ( I  ,  J  )-THET A ( I.J-l)  )* 

1  YBI/DtN 

DFXV=  (  (  V  (IX.P-l)-V  ( I , J-1))*YA-(V  C  I  *  J  ) -V  (I , J-l ) ) »YE ) /CEN 
OF  XS  =  C  (  s  (IX.P-D-S  ( I ,J-1) )*YA-(S  (I,J)-S  (  I  ♦  J-l)  )*YB)/CEN 
D  F  XH -  ( (AH(IX.P-l) -AH { I , J- 1  )  )  *Y A- { Ah ( I  »  J  )  -AH ( I , J-  1 ) ) • YE ) /C EN 
DFYT=(XA«(THETA( I , J)-THETA( I ,J-1) )-( THETA ( I  X , P- 1  ) -THETA ( I.J-l) )*XB 
1  ) /DEN 

DFY  Y=  (  XA«  (  V  (  I  *  J  )  —  V  n,J-l))-(V  C  I  X  -  P-  1 )  -  V  { I , J-l )  )*XE )/CEN 
DFYS=(XA#(S  (I.J)-S  ( I » J-l ) ) - ( S  (IX.P-l)-S  ( I » J-l )  )*XE )/CEN 
OF YH= (XA*( AH( I, J ) - AH ( I  *  J  —  1 )  )  —  ( A  H ( I  X  »  P- 1 ) - AH  *  I  *  J  —  1  )  )*X6)/CEN 
THETA ( I , J )=THETA( I , J-l ) +CFXT*XC+CFYT*YC 
V  (I,J)=V  ( I , J-l ) +CFXV*XC+CFYV*YC 

S  (I,J)=S  ( I ,J-1 )+CFXS*XC+CFYS*YC 

AH  (  I  »  J  )  =  AH  ( I , J-l )+CFXH*XC+CFYH*YC 

I F ( V (  I  ♦ J  )  .GT . 1 . )  GO  TO  440 

AMI  I ,  J)=S0«T(2./(GAM-1.  )  )  *  SORT  ( V  ( I  ,  J  )  **2/ I).  .-V  1 1  ,  J  )  **2 )  ) 

28'.  Y(I,J)=YPM 

X ( I  ,  J  )  =XP ( I ,J) 

PM  (  I  . .! )  =  (  Y  (  I  *  J  )  —  Y  (  I »  J—  D  )  /  (  X  (  I  .  J  )  -  X  (  I. J-l)  ) 

GP=GAM+1. 

GM= (GAM-1, ) /?• 

A  M  2  =  A  M  (  IX.P)**2 
DL-0 . 

I F ( AM ( I ,J)**2-1..LT.0«  )  AM( I  ,J)  =  1.2 
290  ANS  =  GM/(GAM*AM2)+GP/(2.*GAM*AM2)*((1.4GM*AM2)/ 

1  (  l.  +  GM*AM(I  ,  J ) **2 ) )**( G AM/ ( GAM-1. ) ) * ( 1 .- ( GAM*AM  (  I,J)**2)/ 

1  SQKTIAMI I , J)*«?-l. )*DL+DL**2*GAM*AM( I , J ) **2* ( GP*AM ( I , J ) **4-4 .* ( AM 
1  C I • J ) **2-1. ) )/( (AMI  I ,J)**2-1.)**2*4. ) )*EXP(S(  IX,P)-S( I.J) ) 

IF ( ANS.LT.O. )  GO  TO  430 

EPI ( I.J)=ARSIN(SQRT( ANS) ) +THETA ( IX , P ) 

3J  ABC=ABS( (S( I,J)/S(IX,P)-1.)*100.) 

305  I F ( A8C.LT. YT)  GO  TO  410 
GO  TC  370 

31.  GO  TO  (320,330.  KA 

32u  IF(X(I,J).GT.ALP)  GO  TO  36C 

33  CONTINUE 

I F ( (NU( IX)-1).GT.N)  GO  TO  47C 
DO  340 1 J» 1* II 
340  XX ( I ,  I  J  )  =  10.E+10 
NP  T  =  0 

IF(MC.GT.O)  CALL  SHOCK 
IF(NPT.EG.l)  GC  TC  10 
IF(M.GT.C)  GC  TC  35C 
IF(MN.GT.O)  GO  TO  350 
CALL  CSHCCK 
35;  PETUKN 
36c  K A=2 
37  j  N  =  J 
38'  M  =  1 

M AXN=N 
L I  =  I  I-I 
DO  3 8d I J  =  I , L I 
385  XX (  I  ,  I  J  )  =  10.E+ 10 
GC  TO  35*1 

393  AMI  I , J)  =  ( AM(  I - 1 , J)4AM( I - 1 , J  +  l) )/2. 
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400 


* 


4lJ 


420 

430 

440 

450 


460 


470 


480 

490 

500 

510 

520 

530 


GO  TO  200 
N  =  J- 1 
GO  TO  38C 
MC  =0 

WRI T E ( 6  » 500  )  ABC 
GO-1  • 

GO  TO  310 

IF(IR.LE.KK)  GO  TO  400 
WR I TE ( 6 , 5  30 )  IR,J 
GO  TO  400 

EPI  ( I » J )  =  THETA ( I-ltJ  +  1) 

GO  TO  3C0 

AM! I ,J)=( AM( IX,P-1)+AM( IX,P))/2. 

GO  TO  2 80 
UX*U ( AM ( I t J-l )  ) 

UW*U( AM( IX, P)  ) 

AX=A!UX,V ( I » J-l)  ) 

A W= A ( UW , V ( I  X  ,  P  )  ) 

BX=B( JAM,UX, THETA ( I, J-l) ,Y< I ,J-1) ) 

ANX*AN!UX, THETA! I, J-l) ) 

CW=C ( UX » GAM  » RS  ) 

SMAX=SMA!UW, THETA ( IX,P) ) 

CX*C!UW,GAM,RS) 

DW=D (JAM,UW»THETA( I X,  P ) »Y( IX  ,P )  ) 

ZW*Z ( UW t THETA (  IX, P) ) 

TX*T (UX,THETA( I,J-1J ) 

1>X*T  AN  (  THET  A  (  I  ,J-1)+UX) 

HX*TAN(THETA( IX,P)-UW) 

GO  TO  180 

AZ*AX 

AY  =  AW 

BZ  =  BX 

CY=CW 

CZ=CX 

DY=DW 

zz=zw 

TZ  =  TX 
GZ  =  GX 
HZ=HX 
GO  TO  230 
LIM=LIM+1 
J*l  I M 
N=N  + 1 
MAXN*N 
GO  TO  10 

FORMAT (10X,4HYPM( ,13, 1H, 13, 1H) ,E20.8) 

FORMAT! 10X,4HXPM( ,I3,1H, I3,1H) ,E20.8) 

FORMAT (20X,E18,8) 

FORMAT ( 5X ,4  1 10 ) 

FORMAT  ( 1 X » 1 3E 10, 4  ) 

FORMAT! IX, 29HWE  HAVE  NOT  CONVERGEC  FOR  ROW , I  5 , 5HP0 INT  ,  I  5  ) 
END 
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SUBROUTINE  C jHCCK 

COMMON  v'  (  5 , 7  o  )  ,THFTA<5,  70)  ,AM<  5,70  ),S  (5,70  ),X(  5,70  )  ,Y(  5,70  ,NU  5) 
1  ,  E  P I  (  5  ,  )  ,  AH  (  5 , 70  ) ,  PM  (  5  ,  7  j  ) ,  LOC  C  5  )  ,  XP  (  5 , 7C )  ,  X X  (  c  ,  1  5C  )  , 

1  GAM,  RS  » AM  I »  N ,  I  ,  J  ,  M  ,  J  AM ,  ALP ,  TBET  ,  NP  ,  K  A  ,  M  AXN  ,  MN  ,  MC  ,  I M  ,  N  P  T  ,  LF  SH  , 

I  LL.TL.Pl, EMI ,YT, ICEL, IREC.KK.AL, A4,e4,C4,C4,E4,F4, IR 
J  =  N 

299  UX=U( AMI  I , J-l)  ) 

AX  =  A  I UX  »  V l I, J-l)  ) 

BX  =  BI JAM.UX, THETA  I I , J-l  )  ,Y(  I , J-l  )  ) 

CX  =  C IUX.GAM.RS  ) 

GX  =  TAN I THET  A (  I  ,  J-D+IX  ) 

FX  =  T  AN  I EP I (  1-1 , J-l  )  ) 

GM1=GAM-1 . 

AM  I S- AM  I ««2 
GM=GM1/2.*AM IS 
GP1=GAM+1. 

BET=EPI( 1-1 , J-l ) -THETA ( I - 1 , J-l ) 

SA  =  -V(I-1,J-1)*SIMBET)*CCS(EET)*(FX/TAN(6ET)  +  TAMBE7)/FX- 
1  ( 2 . *GM1 ) /GP 1 5 

CGM=(  AMI#SIMEPI  (  I-1*J-1)))««2 
C0MS=(CCM-1, )  *  *2 

SB=-SIN12.*THETAI I-l.J-l) )/SIN{2.*EPI ( I - 1  * J-l) )  +  ( ( G P 1 • A r 1**2 
1  *COM*S INI  THETA (  I-l.J-l)  )**?)/COMS) 

SC=RS/FX«CGMS/ ( ICCM-(GM1/(2»*GAM) ) ) *  I  COM* l GM 1 /2 . )+l. ) ) 

L  =  0 

31 C  L=L+ 1 

FG=FX-GX 

X(  I,  J)  =  ( Y(  I , J-l )-Y( I -  1 , J-l )+FX*X( 1-1, J-l )-GX  *X( I , J-l) )/FG 

Y(  I  ,  J)=(FX*GX*IXI  1-1  ,  J-l) -XU,  J-l  )  )+FX*YC  I»J-1)-GX*Y(  I-  1 ,  J-l  )  )/FG 

DELX*X( I , J)-X( I , J-l) 

EPI C I, J ) =  E  P I  I  1-1 , J-l )  +  ( AX*( V  C I , J -  1 1 - V ( I - 1 , J-l) )-TFETA( I , J-l )  + 

1  THETA ( I-l.J-l )+CX*ISt I , J-l) -S (1-1, J-l) )+BX*DELX)/(SA#AX-Se+SC*CX) 
SEP-SIMEPI  (  I  ,  J  )  ) 

EM=AMIS*SEP**2 

THETA I  I , J ) = AT  AN  I |AMIS*SINl2.*EPI(I,J))-2.*CCSlEPIII,J))/SEF)/ 
1(2.+AMIS*( GAM  +  CCS(2.*EPI ( I ,J) )  )  )  > 

V  T  =  1 • - ( (4.* (EM-1. )«(GAM*EM+1. ) ) )/(GPl**2*AMIS*EM) 

IFIVT.LT. 0. )  GC  TC  3C0 
VC I,J)=SOKTCGM/C l.+GM) 

I F ( (2.*GAM*EM).LT.GM1)  GO  TC  300 
S( I , J )=RS/GM1*AL0G( ( 2 •*GAM*EM-GMi ) /G P 1 ) -G AM*R S/GM 1* 
lAL0G(GPl*EM/(GMl*EM+2. ) ) 

AH ( I ,  J ) = AH  I I - 1 .  J-l) 

I  F ( V ( I, J) .GT.l. )  GC  TO  3GC 
AMI  I ,J)*SGRT(?./GM1)«S0RT(V( I,J)*«2/( l.-VI I,J)««2) ) 

IFIL.EC.l)  GC  TO  340 

IFISI I , J ) .LE.J.O.CR.Xl I, JI.LE.O.C)  GC  TC  3C0 
IFIAMI I, J l.LT.l. )  GO  Tu  3C" 

320  IF  (  Y(  I  ,  J)  .LT.~>.  .OR.  Y  (  I  ,  J  )  .LT.T6ET  )  GC  TO  3C0 
PMC I,J)=(Y( I,J?-YC I,J-1) )/(X( I,J)-X(  ) 

IFCXI  I , J ) .GT.ALP)  GC  TO  33C 
RETURN 
330  M, :  1 

MAXN  =  \ 

RETURN 
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I 


* 


300  N*  J-l 

GO  TO  330 
340  AZ®  AX 
BZ®BX 

cz=cx 

GZ®GX 
F  Z  =  F  X 

350  UX=U ( AM ( I » J )  ) 

AX®.5*( AZ+A(IX,V( I ,J) )  ) 

BX*.5*(BZ+B( JAM,UX,THETA( I, J),VC k,J> ) ) 
CX*.5*(CZ*C(lX,GA^fRS) ) 

‘  FX-i 5*(FZ  +  TAMEPI ( I , J)  )  ) 

GX=.5*(GZ+TAN( THETA( I  J)+UX)) 

“GO  TO  310 
END 


SUflROUl  INE  CIV  ST 

COMMON  V(  5,70)  .THETA  (5,72)  ,AM  (  5,70  ,S(5,70),XC5 ,70  ,  Y(5,7C  ),Kl  (5) 
1  ,EPI ( 5 , 7C ) , AH ( 5 , 7 “ ) , PM ( 5 , 70 ) , LOC ( 5 ) ,  XP { 5 , 7C ) , X X ( 5 , 1 5C  )  , 

1  GAM*  RS  »  AM  I »  N  ,  I  ,  J  ,  M  ,  J  AM  ,  AL P  ,  T E E  T  ,  N P  ,  K  A  ,  M  AXN  , MN  ,MC  ,  I M  ♦  N  PT  t  L  F  SH  ♦ 

1  ll,TL,Pl,EMl ,YT,  IDEL, IkEC,KK,AL,A4,R4,C4,D4,E4,F4, IR 
J  =  1 

U  X  =  U  ( A  M  (  1-1,2) ) 

AX=A(UX,V(  1-1,?)  ) 

CX-CIUXtGAMfPS ) 

DX=D( JAP , UX, THETA ( 1-1,2) , V ( 1-1,2) ) 

HX  =  TAN(THETA{  I  —  1  *2  )  —  IjX) 

EX=TAN( THETAC  1-1,1)) 

N=N-NP+MC 
M=  1 
NP  =  2 

GM--GAM-1. 

L  =  0 

110  L=L+ 1 

S(I,1)=S(I-1,1) 

AH ( I , 1 ) *  AH ( 1-1,1) 

DEN=EX-HX 

xnxc  i  ,n 

X(I,1)*(Y(I-1,2)-Y(I-1,1)+EX«X(I-1, 1 )-HX«X( 1-1,2) )/CEN 
DELX  =  X ( 1 , 1 )  —  X (  1-1,2) 

Y(  1 ,1)*(EX*HX*(X(  I  — 1,1)  —  X(  1-1,2)  )  +  EX« Y (  I  —  1 ,2  )— f-X*YI  1-1,1)  )  /  DEN 
XQ  =  X ( I , 1 ) 

AM  Cl ,  1)=SQRT(2./GM*{ { EM  1**2 *GM/ 2 . ♦  1 .  )  * ( ( P 1/P2 ( XC , A4 , BA ,C4 , C A , E4 ,F4 
1  )  )*»(GM/GAM)  )-!,)) 

VC  I » 1 ) =SQRT ( AM (  I,l)»«2/(2./GM+AM(I,l)**2)) 

T1=THETA( 1,1) 

THETA ( I , 1 )  =  THET AC  1-1,2 )-AX»(V( 1 ,1)*SCRTCAH( 1,1) /AHC I- 1 , 2 ) ) -V ( I- 1 , 
1  2))+DX«DELX-CX*($(I-l*l)-S( 1-1,2) )+AX«(AH(I,l)-AH( 1-1,2) )/ 

1  CVC 1-1 , 2 ) • <  AH  f  1,1  )  +AH ( 1-1,2) )) 

EPIC I,l)=° 

IFCA6SC (Xl-XC I  ,I))/X(1 ,1 ) ) .IT. .OCOl. ANC.AfiSC  (Tl- THETAC 1 , 1) )/ 

1  THETAC  1,1)  J.LT..CC I'l)  GO  TO  120 
IFCL.GT.5Q)  GO  TO  115 
IF(L.NE.l)  GC  TO  1«<" 

AZ  =  AX 
CZ=CX 
DZ  =  DX 
HZ=HX 
EZ  =  E  X 

180  UX=U ( AM ( I , 1 ) ) 

AX=.5«( AZ+ACLX.VC 1,1)) ) 

CX=.5*(CZ+C(LX,GAM,RS))  -  - 

DX=. 5« (CZ+O (JAM ,UX, THETAC I , 1 ) , Y ( I , 1 ) ) ) 

HX  =  .5#(HZ  +  TAMTHETA(  I,l)-tX)) 
cX=.5*(E7+TAN( THETAC 1,1) ) ) 

GO  TO  llu 
115  WRITE (6,160) 

IZu  IF(Y( I,  1)  .LT.O..OR.Y ( I, 1 ) .LT.TBET)  GC  TC  130 
IF(N.EQ.l)  GO  TC  130 
RETURN 

130  WRITE  (6,  140  I 

WRirE(6,15i)V(I,J),THETA(I,J),AM(IfJ),S(I,J),X(I,J),Y(I,J), 
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1  AH  (I,J) 

N  =  1 

RETURN 

14j  F0RMAT{///3H  ROW  ,I2*15H  TEST  STOP///) 

160  FORMAT ( 7E 18. 8 ) 

160  FORMAT (40X»40HWE  HAVE  NOT  CONVERGED  ON  THE  FIRST  POINT) 
END 


APPENDIX  B 


PROGRAM  FOR  VISCOUS  REGION 

Once  the  subsonic  portion  of  the  boundary  layer  is  subdivided 
into  n  strips,  the  flow  in  each  strip  is  assumed  to  be  governed  by  the  one¬ 
dimensional  equations  with  friction  and  heat  transfer.  The  resulting  system 
of  equations  consist  of  "3n"  nonlinear,  ordinary  differential  equations. 

These  equations  are  connected  through  the  boundary  conditions  on  heat 
transfer,  q,  and  friction,  c^  (Fig.  6).  Once  the  starting  conditions, 
pressure  distribution,  and  shape  of  basic  streamline  are  known,  the 
equations  are  solved  simultaneously  by  the  Runge-Kutta  method. 

At  every  station  station  x  it  is  therefore  possible  to  find  an  area 
consistent  with  the  flow  conditions  in  each  stream  tube  in  the  flow  field. 

Since  the  basic  streamline  and  the  area  for  the  viscous  layer  are  now 
known,  the  location  of  the  Dividing  Streamline  can  be  obtained.  The  point 
at  which  the  Dividing  Streamline  intersects  the  axis  is  the  location  of  the 
rear  stagnation  point. 

Downstream  of  the  rear  stagnation  region  the  inviscid  characteristic 
and  viscous  layer  analyses  must  be  such  that  the  basic  strearruine  is  located 
at  radial  location  where  the  outer  streamline  of  the  viscous  layer  exacdy 
coincides  with  the  inner  streamline  of  the  inviscid  flow  field.  If  at  any 
x  station  this  condition  is  not  satisfied,  then  the  assumed  pressure  along 
the  basic  streamline  is  changed  until  the  above  condition  is  satisfied. 

For  details  of  the  viscous  layer  program,  one  may  consult  the  flow 
chart  (Fig.  B-l)  and  the  program. 
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SUBROUTINE  VISCOUS 


* 


FiG.(B-l)  VISCOUS  PROGRAM 
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COMPUTER  PROGRAM  FOR  SHEAR  LAYER 


SUBROUTINE  RK ( X ) 

DIMENSION  S ( 20 ) ,F(2C) , H ( 20 ) , P ( 20 ) ,U ( 20 ) ,A(50),B( 10  )) , R ( 2 2 ) , SH ( 21 ) 

1  ,DN(20) ,0(21) ,RX(20) , XR ( 20 ) , TX ( 20 ) , AH ( 80 ) , AS ( 80 ) , AF ( 80 )  , C ( 20 )  , 

2  CF120) »FF(20) ,HF(20) ,V( 20) ,CCF(20) 

DHX1(H,F,C) =-DX#( (  AME*SQRT(H)*6.2831853)*PE/(SQRT(F)*PF*C)MGHF* 

1  SF-QLF*RF  )  ) 

DCX(H,F,C )=DX* ( C  * ( ( l.-F )*UF/(GAM#PF*F )+(HX/H)*( l.+GMD*F ) 

1  ♦( ( l.+GM*F)*CFT  /ABS(SF-RF) )  )  ) 

DFX(H,F,C)=DX*(-F*( l.+GMD*F ) * ( 2 . *UF / ( GAM*PF*F ) +HX/H+CF T*? . / A6 S ( SF 
1  -RF)  )  ) 

READ! 5, 150)  NO 

READ! 5, 170)  AME ,REL , ALP , GAM , P E , TH 

READ! 5,155)  DX,AXX,E1,E2 
RE AD (5, 155)  VL,SHO,VO 
Z*0. 

NQD=NQ/3 
READ! 5,170) 

READ! 5, 170) 

READI5, 170) 

P.EAD(5,150) 

1 F ( IY.NE.O) 

WR I TE ( 6 , 270 ) 

READ! 5, 150)  SAME 
LIM=5*NQD 

I F ( SAME*  NE»  0 • )  GO  TO  7 
L  I  M=5 


11 


111 


(C ( I ) , 1=1, NOD) 
(F( I ), 1=1, NOD) 
( H ( I ) , 1=1, NOD) 
IY 

GO  TO  147 


NX  =  1 

READI5, 170) 
LI  =2*L  IM 
READ (5, 170) 
READ ( 5,150) 
READ! 5, 170) 


( A ( I )  ,1  =  1, LIM) 


) 


( B ( I  )  ,  1*1 ,L I 
I  XR 

(XR(I),RX(I),TX(I) , I =1 , I XR ) 

WRI TE ( 6 , 250 )  E1.E2 
WRITE(6,255)A( 1),B(1),B(2),B(2) 

WRITE (6, 260)  ( A ( I -1 ) , A ( I ) , B ( 2* I - i ) , B ( 2* I ) , B ( 2* I ) , I = 2 , L I M ) 
WR I TE ( 6 , 265 )  (XR(I),RX(I),TX(I),I=1,IXR) 

NQDP=NQD+ 1 
NQDM=NQD- 1 
NW=NQD*4 
NX=NQD*5 
PH=3. 14159265 

PCF=2.*PE*C0S( ALP)/SIN( ALP)/REL 
S ( NOD )  =  1  • 

CCFINQD)*0. 

GM=GAM-1 . 

GM0=GM/2. 

HN=1.+GMD*AME**2 
EL=2./(REL*SIN( ALP)*(HN) ) 

XM=-1 . 

I  XF  =  0 

AB=COS (  TH) 

DO  111  1=1, I XR 
IF(X.LE.XR(  I  )  )  GO  TO  112 
C0‘  TINUC 
1  =  1X0 
GO  TO  113 


71 


112  IFIX.EQ.XRtf ).0R.X.LT.XR<1))  GO  TO  113 
RIl)*RXtI-l)+(RXII)-RXtI-l))/tXRtI)-XRII-l))*tX-XR(I-l) ) 

TTH«TXII-1)-MTX!I)-TXII-1))/IXR1I)-XRII-1))*IX-XRII-1)) 

R  C 1 1 *R (II /COS( TTH) 

GO  TO  114 

113  RI1)*RX<I)/C0S<TXII)) 

TTH*TX I I ) 

114  IFIXM.GT.O.)  R(1)»0. 

I F ( IY.NE.O)  R(1)«0. 

XT«X 

IW*1 

DO  12  I»1,NQD 
CFI I )«C( I ) 

FFII  )*FU  ) 

12  HF 1 1 ) *H( I ) 

ZF*Z 

IF (XM.LT.O. )  GO  TO  141 
IFIX.EQ.O.)  Z«0. 

25  SHINQDP )*E1+E2*Z 

00  70  I *1 »NQD 
IF (HFl I ) .LT.O. )  RETURN 
IF (FF II ) • LT.O. )  RETURN 
SMI  I  )*HN*HFI  I  )/l  l.*GMD*FFU  )) 

VU)*SQRT(FF(I)«SH(I  H/AME 
DNII)«  SIII-RII) 

IF ( IY.NE.O)  SHlNQOPIxSHO 

70  Q(  I  )»EL*SQRT(SH(  1-1)  )*ISHI I  J-SMI  1-1)  )/(DN(  I-U+DNI I  ) ) 

Q ( NQDP ) «EL*  SORT (SHI NQO ) )*I SHINQDP )-SH  I NQO) ) / (ONI NQO ) *2. ) 

I F ( XM.GT.O. )  Q(NQ0P)*O. 

IF ( IY.NE.O)  Q(NQDP)*0.0 
I F ( IY.NE.O)  QI  I  )**-Q(  I ) 

QIl)»EL*SQRT(SH0)*(SH(l)-$H0)/tDNli)*2.0) 

I F ( IY.NE.O)  Q(l)-0.0 
IFIIXF.LT. 1)  GO  TO  1391 
IFIIH.EQ.3)  GO  TO  120 
K«1 

Z«ZF+DX#COS(TTH) 

L*0 

00  100  I*l*NXf.5 
L*L  +  i 

I F ( Z.LE.A ( I ) )  J*K 

IFIAII  ) .LT.Z.AND.Z.LE.AI I  +  l)  )  J*K«2 
IF(A(I+1).LT .Z.AND.Z.LE.At 142) )  J*K«4 
IF (A(  1+2)  .LT.Z.AND.Z.LE.AI  1+3) )  J«M!6 
IF (A( 1+3) .LT.Z.AND.Z.LE.AI I+A) )  J*K+8 
P(L)»B(J)+B(J+1)*Z 
95  U( L) *B( J+l ) 

100  K-K+10 

IF! SAME.NE.O. )  GO  TO  120 
00  110  I«2.NQD 

pi  n-pci) 

110  U(I)«U<1) 

IF! IY.NE.O)  VL*VO 
120  00  125  I *1 » NOOM 

123  CCF(I)»PCF*SH(I)**l.5*(V(i  W(H-*))/(P(I  )*V(I  )**2*.5*!R(I)-R(I+2) 
1  )) 
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CCF{NQD)=PCF*SH(NQD)**1.5*(V(NQD)-VL)/(P(NQD)*V(NQD)* 

1  { ft ( NOD ) -R ( NOOP ) ) ) 

I  =0 

00  130  J  =  I W » NW , 4 
1  =  1  +  1 

CCFO«PCF*  SHO#* 1.5*1 VO-V ( 1 ) ) / ( P ( 1 ) *V0**2* ( R ( 1 »— R 12)  ) ) 

CFT  =  +  CCF  (  I )  —  CC F  (  1-1). 

IF(I.EQ.l)  CFT=+CCF( 1)-CCF0 
IFUY.NE.O)  CCF0=0.0 
IFUY.NE.O)  CF  T =-CFT 
PF  =  P {  I ) 

QHF*Q ( 1+1 ) 

QLF-Q ( I ) 

RF  =  R  (  I ) 

SF-SII  ) 

UF=U (I  ) 

IF(FF(I).LT.O.)  RETURN 

AH ( J ) *OHX 1 ( HF ( I ) ,FFt I ),CF(  I )  ) 

HX=AH( J)/DX 

AS(J)  =  OCX(HFU),FF(I),CF(I)) 

AF(J)*DFX(HF( I )»FF< I ) ,CF(I  )  ) 

130  CONTINUE 

GO  TO  (132* 135* 136 • 138 ) • IW 

132  I W=2 
X*XT+DX/2  • 

IFIXM.GT.O.)  R(1)«0. 

J«i 

0  =  2. 

133  DO  134  1=1, NOO 
CF(I)*AS(J)/0  +C(I) 

FF( I )*AF(J)/0  +F( I  ) 

HFU)=AH(J)/0  +H(I) 

TC—CFII  ) 

IFUY.NE.O)  TC*CF  (  I  ) 

I F ( ( TC/PH+R (I)**2).LT.j.)  RETURN 
S( I ) =SQRT ( TC/PH  +R{ I )**2) 

R ( 1  +  1 ) =S(  I ) 

134  J= J+4 

GO  TO  25 

135  I W*3 
J  =  2 

GO  TO  133 

136  I  W=4 
X-XT+DX 

IFIXM.GT.O.)  R ( 1 ) =0. 

J  =  3 
0  =  1. 

GO  TO  133 
138  J  =  1 

DO  139  1=1, NQO 

C(I  )=C(  I)+<AS(J)+2.*AS(J+l)+2.*AS(J+2)+AS(J+3)  )/6. 

H  (  I  )  =H  (  I  )+{AH{  J  )+2.*AH<  J  +  l)  +  2.*AH(.»  +  2)  +  AHC  J  +  3)  )/6. 
FII)=F(  I  )  +  (AFU)+2.*AF(J+l)  +  2.*AF(,H2)  +  AFI  J  +  3)  )/6. 
TC=-C(I) 

IFUY.NE.O)  TC=C  (  I  ) 

RT=TC/PH+R< I )**2 
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1385 


139 

1391 


140 


141 


142 

143 


144 

147 

150 

155 

170 

180 

210 

250 

255 

260 

265 

270 

275 


IFIRT.GE.O. )  GO  TO  1385 
X=500. 

GO  TO  139 
S ( I ) =SQRT ( RT ) 

R { 1+1 ) =S ( I ) 

J= J+4 
I  XF*2 

X1=X»C0S !  TTH ) 

WRITE!6,210)  X,Z 
00  140  1  =  1 » NQOP 
IF(F(I».LT.O.)  RETURN 
FS=SQRT!F( I ) ) 

WR I TE ( 6 » 180 )  C(I) »H(I) ,SH(  I  ) ,Q(  I  )  ,FS  ,CCF{I),R(I) 

CONTINUE 

IF(Z.LT.AXX)  GO  TO  11 
RETURN 

00  142  1*1, NOD 
TC*-C ( I ) 

IF ( I Y.NE.O)  TC=C (  I  ) 

ZT=TC/PH+R! I )**2 
IFIZT.LT.O.)  GO  TO  143 
SI I)*SQRT!ZT) 

RtI*l)»S( I) 

GO  TO  25 
XM=1  • 

AB=1. 

DO  144  1=1, NOD 
S< I)=Sl I ) *COS ( TTH ) 

R(I  +  l)*Sm 

C(I)=PH*(SII)**2-R(I)**2) 

GO  TO  11 
WR I TE ! 6 , 275 ) 

GO  TO  5 
FORMAT! I?) 

FORMAT ( 5E 15*  8 ) 

FORMAT ( 6E 13*  6 ) 

F0RMAT(3X,7E18.8) 

FORMAT! ///, 2X,2HX=,2E15.8,/12X,1HC,17X,1HH,16X,2HSH,17X,1FQ,16X, 

1  3HSQF, 16X,2HCF,16X,1HR) 

FORMAT!  10Xt2HH=,E15.8,lH+,E15.8,2H»X,//) 

FORMAT (26X,8HIF  Z  LE  ,E13.6,5H  P= , E13.6, 1H+, E13.6V 9H*Z  AND  t=y 
1  E13.6) 

F0RMAT(5X,2HIF,E13.6,14HLT  Z  AND  Z  LE  ,E13.6,5H  P* , El 3. 6 , 1H* , 

1  E13.6,9H*Z  AND  U=,E13.6) 

F0RMAT(///,21X,1HX,19X,1HR,17X,5HTHETA,/, !10X,3E20.8) ) 
F0RMAT(1H1,25X,72HW6  ARE  WORKING  FROM  TRAILING  EDGE  OF  CCNE  ( X=iO ) 
1  TO  REAR  STAGNATION  POINT,///) 

FORMAT! 1H1,35X,  56HWE  ARE  WORKING  DOWNSTREAM  OF  REAR  STAGNATION 
1  POINT  ( X=0) ,///  ) 

END 
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